“ — 

/  AD-Allb  774 

UNCLASSIFIED 

AIR  FORCE  INST  OF  TECH  WR2 SHT-PATTERSON  AFB  OH  F/6 

LABYRINTH  SEAL  EFFECTS  ON  ROTOR  BEARING  SYSTEM  STABILITY. (U) 

MAY  62  A  U  PAVELKO 

AFIT/NR/62-10T  N, 

3/9 

ii 

tm 

1 

III 

II 

II 

i 

HI 

i 

HI 

|| 

' 

i 

1  M  \i 

mm 

MB 

■L 

■ _ 

Dre  FILE  COPY  _  ad  All 6774 


nisiri  ass 


© 


REPORT  DOCUMENTATION  PAGE 

READ  INSTRUCTIONS 

BEFORE  COMPLETING  FORM 

1.  REPORT  NUMBER  J  2.  GOVT  ACCESSION  NO. 

AF IT/NR/82- 10T  /  /\£)_  fr  /  j  ^  99 

r . . 

4.  TITLE  (and  Subtitle) 

Labyrinth  Seal  Effects  On  Rotor  Bearing  System 
Stability 

5.  TYPE  OF  REPORT  &  PERIOD  COVERED 

THESIS/D mtmw 

6.  PERFORMING  ORG.  REPORT  NUMBER 

7.  AUTHORfjJ 

Anthony  John  Pavel ko 

8.  CONTRACT  OR  GRANT  NUMBERfsJ 

».  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

AFIT  STUDENT  AT:  University  of  Texas  at  Austin 

10.  PROGRAM  ELEMENT.  PROJECT.  TASK 
AREA  &  WORK  UNIT  NUMBERS 

1.  CONTROLLING  OFFICE  NAME  AND  ADORESS 

AF IT/NR 

4PAFB  OH  45433 

12.  REPORT  OATE 

May  1982 

13.  NUMBER  OF  PAGES 

170 

14.  MONITORING  AGENCY  NAME  &  AODRESSf//  different  from  Controlling  Office) 

15.  SECURITY  CLASS,  (of  this  report) 

UNCLASS 

15*.  DECLASSIFICATION/ DOWN  GRADING 

SCHEDULE  | 

_ _ _ f  r 

APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


< 


17.  DISTRIBUTION  STATEMENT  (of  the  abstract  entered  In  Block  20,  H  different  from  Report ) 


N, 


is.  supplementary  notes 

APPROVED  FOR  PUBLIC  RELEASE:  IAW  AFR  190-17 


r/# — . 

^CpfN  E.  WOLAVER 


22  JUN  'aS2 


19.  KEY  WOROS  (Continue  on  reverse  side  if  necessary  and  Identity  by  block  number) 


Dean  for  Research  and 

Professional  Development 
AFIT.  Wriqht-Patterson  AFB 


50.  ABSTRACT  (Continue  on  reverse  side  If  necessary  and  Identity  by  block  number ) 


ATTACHED 


"07  07  °59 


OH 


DD  I  JAN  73  1473  COITION  Of  I  NOV  «S  IS  OBSOLETE 


UNCLASS 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (When  Dele  Entered) 


ABSTRACT 


This  thesis  presents  a  stability  analysis  of  rotor  bearing 
system  operation  affected  by  a  single  labyrinth  seal.  Labyrinth  seals, 
primarily  intended  for  leakage  control  in  turbines,  have  been  observed  to 
affect  rotor  stability  in  high  pressure  steam  turbines.  A  small 
experimental  rotor  bearing  system,  existing  in  the  Mechanical  Engineering 
Department,  and  a  compatible  labyrinth  seal  are  the  subjects  of  this 
analysis.  Linear  analytical  models  of  these  rotor  bearing  and  labyrinth 
seal  systems,  both  isolated  and  coupled,  are  developed.  Labyrinth  seal 
geometry  and  pressure  variations  generate  various  affects  on  rotor 
bearing  system  stability.  These  affects  can  help  develop  labyrinth  seal 
design  criteria  directed  toward  high  pressure  turbine  stability 
enhancement . - 
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NOMENCLATURE 


a  Seal  strip  dimension  (in) 

2 

A  Area  (in  ),  Matrix  descriptor,  Variable,  Seal  position  identifier 
b  Damping  coefficient  (lfb  sec/in),  Seal  strip  dimension  (in) 

B  Damping  coefficient  (lbf  sec/in).  Variable,  Seal  position  identifier 

c  Bearing  to  journal  radial  clearance  (in) 

C  Seal  position  identifier,  Variable 
^  Function 

e  Eccentricity  (in),  Natural  logarithm  base  =  2.718 
E  Modulus  of  elasticity  (lbf/inZ) 
f  Friction  factor  (dimensionless) ,  Function 
F  Force  (lbf) 

h  Minimum  bearing  to  journal  or  seal  base  to  disc  radial  clearance  (in) 
H  Related  to  h 

I  Moment  of  inertia  (lbf  in  sec2) 

^  Stffftttisss  factor  (lbf /in).  Specific  heat  ratio  (dimensionless) 

K  Stiffness  coefficient  (lbf /in).  Constant  coefficient 
l  Seal  chamber  length  (in),  Journal  length  (in) 

Lspt.  Rotor-bearing  system  overall  length  (in) 
m  Mass  (lbf  sec2/in) 

2 

M  Mass  coefficient  (lbf  sec  /in),  Mach  number  (dimensionless) 

n  Exponent  in  seal  friction  relations,  Integer 

0  Fixed  reference  point 

2 

P  Pressure  (lbf/in  ),  Momentum  (lbf  sec) 
q  Generalized  coordinate 


viii 


r  Disc  radius  (in).  Journal  radius  (in) 

R  Gas  constant  (in^/sec^°R) 

Re  Reynolds  number  (dimensionless) 

s  Seal  circumferential  distance  independent  variable  (in) 
t  Time  (sec) 

T  Temperature  (°R) ,  Kinetic  energy  (lbf  in) 
u  Seal  circumferential  fluid  velocity  (in/sec) 

U  Disc  or  journal  surface  velocity  (in/sec) 
v  Eigenvector  element  (dimensionless),  Velocity  (in/sec) 

V  Potential  energy  (lbf  in) 

W  Weight  (lbf) 

x  Variable,  Coordinate 
y  Variable,  Coordinate 

Y  Related  to  y 

z  Variable,  Coordinate 
Greek 

A  Finite  increment 

X  Seal  angular  position  variable  (radians) 

S  Variation 

E.  Eccentricity  ratio  (dimensionless)  =  e/c 
©  Bearing  angular  position  variable  (radians) 

^  Eigenvalue  (dimensionless) 

2 

/a,  Dynamic  viscosity  (lbf  sec/in  ) 

^  Generalized  force  (lbf) 

2 

»  Kinematic  viscosity  m  IT/P  (in  /sec) 
if  Constant  5^3. 14 


r 


^  Density  (lbf  sec/in  ) 

2 

T  Frictional  shear  stress  (lbf/in  ) 

4>  Journal  attitude  angle  (radian) 

ti)  Angular  velocity  (rad/sec) 


Subscripts 

letters  Identification 


numbers  Identification 


Seal  location,  Axial 


Seal  location.  Bearing 
Seal  location,  Chamber 


Eccentricity  direction 


Related  to  friction 


High  clearance  area,  Related  to  bearing  or  seal  minimum 
clearance 


Inside,  Particular  value 


General  value 


Journal 


Low  clearance  area,  Left 


Nominal  value 


Outside,  Initial  value 

Related  to  seal  chamber  pressure 

Right 

Seal,  Steady 

Stiffness,  Seal  circumferential  distance 


X 


T  Total,  Differentiating  subscript 

u  Related  to  seal  circumferential  fluid  velocity 

W  Whirl,  Discharge  coefficient 

x  Coordinate  frame  direction 

y  Coordinate  frame  direction 

2  Coordinate  frame  direction 

Greek 

^  Related  to  journal  attitude  angle 
Other 

_  Vector  or  matrix 

dt 

d*r 

Other  value,  (  T  =  tU/r) 
Coenergy 
Exponentiation 


Superscripts 

•  • 

/ 

* 

numb er / va r iab 1 e 


Chapter  1 


INTRODUCTION 


Steam  Turbine  (Figure  1. 1)  instability  in  the  form  of 
precession  or  whirl  (  (j>  )  of  a  rotor  about  a  fixed  reference  position 
(Figure  1.2  -  Turborotor,  circumferential  cross  section)  has  received 
extensive  study  throughout  this  century.  Such  instability  is 
characterized  by  cross  coupled  orthogonal  forces  acting  on  a  turborotor 
to  promote  rotor  precession  (Figure  1.2)  M  .  Such  instability  is 
generally  classified  as  either  forced  or  self  excited. 

Forced  rotor  whirl  is  caused  by  mass  unbalance.  It  is 
characterized  by:  (1)  a  whirl  frequency  that  equals  rotor  rotational 
frequency,  (2)  whirl  amplitudes  that  peak  within  a  narrow  speed  band 
(Figure  1.3  -  Forced  vibration  attributes  in  Rotating  Machinery),  and  (3) 
an  absense  of  oscillatory  rotor  fiber  stress  ^2^  .  This  form  of 

instability  is  usually  corrected  through  the  use  of  balance  weights. 
When  balance  weight  application  does  not  restore  turborotor  stability, 
self  excited  vibrations  exist. 

Self  excited  vibrations  generally  are  caused  by  flow  or 
friction  energy  that  generates  rotor  whirl  [3]  .  Such  vibrations  are 
characterized  by:  (1)  a  whirl  frequency  that  is  nearly  constant  and 
independent  of  rotor  rotational  speed  and  is  at  or  near  a  rotor  natural 
frequency,  (2)  whirl  amplitudes  that  suddenly  increase  as  a  particular 
rotor  rotational  speed  value  is  achieved  and  continue  to  increase  with 
increasing  rotor  speed,  (3)  alternating  stresses  in  rotor  fibers,  caused 


by  nonsynchronous  whirl,  that  make  self  excited  whirl  more  destructive 


than  forced  whirl 
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Because  this  instability  form  is  as  complicated  as  it  is  destructive  it 
requires  additional  attention. 

Instability  phenomena  such  as  Oil  Whip  and  steam  excited  whirl 
are  common  forms  of  self  excited  rotor  whirl.  Oil  Whip  is  journal 
bearing  induced  rotor  whirl  that  occurs  when  a  pressurized  bearing  film 
thickens  and  begins  to  trail  the  position  of  minimum  journal  to  bearing 
clearance.  (Figure  1.4  -  Oil  Whip) [4] . 

This  trailing  high  pressure  fluid  region  imposes  a  force,  Fw, 
on  the  journal's  centroid,  normal  to  rotor  eccentricity,  e  (00'  in  Figure 
1.4),  and  in  the  direction  of  rotor  rotation.  Rotor  whirl  at  a  frequency 
of  h  that  of  rotor  rotational  speed  (u>)  results.  This  whirl  is  further 


h  -  minimum  journal  to  bearing  clearance 


Figure  1.4  -  Oil  Whip  Transverse  cross-sectional  view 


defined  as  "forward"  whirl  since  precession  ( ^  ) 
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occurs  in  the  same 
direction  as  rotor  rotation.  Oil  Whip  occurs  when  rotor  rotational 
frequency  is  increased  beyond  a  system  particular  value  known  as  the 
"Stability  Threshold"  .  Although  improvements  in  bearing  design  have 
increased  this  threshold  value  in  many  rotor  bearing  systems,  steam 
turbines  must  still  contend  with  other  sources  of  self  excited  whirl. 

Self  excited  whirl  energy  transformed  from  steam  flow  has 
increasingly  affected  steam  turbine  stability  with  the  advent  of 
increasing  upstream  turbine  pressures.  Whirl  forces  (Figure  1.2), 
necessary  to  cause  "Steam  Whirl",  are  produced  by  unbalanced  torques 
This  torque  unbalance  results  from  higher  turbine  blade  energy  losses  at 
high  shroud  to  blade  row  clearance  areas.  Caused  by  an  initial  rotor 
deflection,  steam  leakage  at  these  large  clearance  areas  will  exceed 
leakage  at  opposite  smaller  clearance  areas  (Figure  1.5).  Consequently, 
the  turbine  blades  nearest  these  high  clearance  annular  sections  will 
receive  less  fluid  energy  than  will  opposite  small  clearance  annular 
section  blading.  The  resulting  unbalanced  torques  produce  a  force  that 
acts  on  the  rotor's  centroidal  axis  and  is  orthogonal  to  rotor 
eccentricity  (e) . 

Nonsynchrous  forward  whirl  at  the  fundamental  natural  frequency 
of  the  rotor  results.  This  whirl  form  is  load  dependent  because  degree 
of  torque  unbalance  increases  with  increasing  upstream  pressure.  Since 
Steam  Whirl  is  a  relatively  new  phenomenon  present  correction  methods  are 
limited  to  bearing  modification  and/or  rotor  stiffening  CO  .  These 


techniques,  such  as  the  use  of  tilting  pad  bearings,  often  introduce  new 
problems.  For  example,  tilting  pad  bearings  provide  a  reduced  defense 
against  mass  unbalance  instability  with  respect  to  standard  cylindrical 
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f^  =  flow  through  high  clearance 
Shroud  areas 


Blade  Row  OD  f  =  flow  through  low  clearance 

areas 


Figure  1.5  -  Steam  Leakage  Between  Turbine  Shroud  and  Blade  Row 

journal  bearings  M  .  Therefore,  alternate  stabilization  methods,  with 
minimal  effects,  are  needed  to  enhance  whirl  protection  for  high  pressure 
steam  turbines. 

One  source  of  such  whirl  protection  has  existed  since  the 
inception  of  the  steam  turbine.  It  is  derived  from  cross  coupled 
orthogonal  forces  generated  by  steam  flow  through  labyrinth  seals. 
Labyrinth  seals  were  originally  intended  to  reduce  steam  leakage  between 
rotor  and  stator  turbine  elements  by  providing  flow  resistance  M  • 
Alford  suggested  that  these  seals  were  a  source  of  instability  when 
rotor  deflection  caused  seal  flow  area  at  a  longitudinal  cross  section  to 
be  converging  (Figure  1.6b).  Correspondingly,  he  also  believed  that 
divergent  seal  flow  area  enhanced  rotor  stability  (Figure  1.6c). 


a) 


Turborotor 
Cross-section  at 
Labyrinth  Seal 


Stator 


Labyrinth 

Seal 


4 


b)  View  A-A 

Convergent  Seal 
(h-a)>  (h-b) 


[e-j Z  -»| 


Flow  Direction 

- > 


c)  View  A-A 

Divergent  Seal 
(h-a)<  (h-b) 


!j  Rotor 

r 

Q 

K 

Stator 

Figure  1.6  -  Labyrinth  Seal  Geometries 
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Unfortunately,  his  analyses  neglected  circumferential  flow  through  the 
seal  cavity.  Spurk  and  Keiper  [•]  made  a  claim  directly  opposite  to 
Alford's  —  that  divergent  seal  geometry  promoted  instability. 
Unfortunately,  their  analysis  was  also  incomplete.  They  neglected  rotor 
rotation.  Because  of  such  disagreement,  as  well  as  potential  for 
stability  enhancement,  labyrinth  seal  forces  require  more  thorough  study. 

Shatof f  [0  provides  an  explanation  supporting  Alford  while 
still  suggesting  that  divergent  and  convergent  seal  geometries  can  be 
either  stabilizing  or  destabilizing  depending  upon  other  operating 
conditions.  Figure  1.7a  depicts  convergent  seal  geometry.  Here,  Fg  and 
F^  are  components  of  the  resultant  force  produc  ed  by  the  seal  (F^  )  • 
Since  the  whirl  component  of  F^s  (F4  )  exceeds  the  bearing  damping  force, 
forward  whirl  at  the  rotor's  fundamental  natural  frequency  results. 
Alternately,  divergent  seal  geometry  (Figure  1.7b)  induces  backward  whirl 
at  an  otherwise  stable  rotor  speed  (a)).* 


This  suggests  that  "negative"  bearing  damping  (reversed  br^ 
vector  direction  with  respect  to  Figure  1.7a),  associated  with  rotor 
speed  increase  beyond  the  threshold  of  Oil  Whip  can  be  offset  (Figure 
1.7c).  Therefore,  by  increasing  the  Oil  Whip  stability  threshold, 
divergent  seal  geometry  enhances  rotor  bearing  system  stability. 
Similarly,  since  Steam  Whirl  is  a  forward  whirl  phenomenon  that  enhances 
oil  whip  instability,  it  too  might  be  counteracted  by  backward  whirl 
inducing  divergent  seal  geometry.  Although  seemingly  plausible, 
experimental  verification  is  needed  to  prove  such  a  theory.  Wright's 
development  and  operation  of  a  seal  and  disc  test  apparatus  partially 


*  As  defined  by  Oil  Whip 


a)  Convergent  Seal 


b)  Divergent  Seal 


c)  Divergent  Seal 
Beyond  Oil  Whip 
Threshold 


Fbs= 


F  = 
e 


* 

br4 


Resultant  Bearing 
Fluid  Force 

k  e  -4-  k,  e 
s  b 

Shaft  and  Bearing 
Elastic  Forces 

Static  (elastic) 
Seal  Force 

Dynamic  (damping) 
Seal  Force 

=  Bearing  Damping 
Force 


Figure  1.7  -  Labyrinth  Seal  Forces 
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satisfies  this  requirement. 

Wright's  model  was  not  representative  of  actual  high  pressure 
steam  turbines.  Its  8-inch  diameter  disc  and  .05  inch  nominal  seal  strip 
to  disc  clearance  were  much  smaller  than  their  high  pressure  steam 
turbine  counterparts.  Additionally,  he  did  not  attempt  to  accurately 
scale  rotor  surface  speed  and  friction  factors.  However,  his 
experimental  efforts  did  provide  a  basis  for  follow-on  experimentation 
with  actual  high  pressure  turbine  seals. 

Wright's  work  also  provides  verification  of  Alford's  theory 
that  divergent  seal  geometry  is  stabilizing.  Therefore,  his  results  are 
probably,  at  least,  qualitatively  reliable  and  warrant  some  discussion. 

With  divergent  seal  geometry,  Wright  produced  strong,  self 
excited  backward  whirl  —  an  effect  possibly  capable  of  counteracting  the 
forward  whirl  of  Steam  Whirl  and  Oil  Whip.  Although  his  model  was  not 
ideally  suited  for  accurate  convergent  seal  geometry  experimentation,  his 
results  did  indicate  reduced  forward  whirl  stability  for  convergent 
geometry  relative  to  divergent  and  straight  (a=b  in  Figure  1.6) 
geometries.  Such  consistency  and  apparent  plausibility  encouraged  him  to 
consider  using  his  results  as  a  basis  for  alternate  turbine  seal  design. 
He  believed  that  if  seal  geometry  affected  rotor  stability,  then  design 
of  labyrinth  seals  to  enhance  stability  was  possible. 

Such  design  information  must  provide  formulas  for  seal 
parameters  and  materials  and  be  applicable  for  the  high  pressures  and 
speeds  of  today's  high  pressure  turbines.  Analytical  modeling  is 
required  to  provide  such  mathematical  relations.  The  first  step  in  such 
modeling  is  analytical  verification  of  Wright's  results,  particularly, 
that  of  stabilizing  divergent  seal  geometry.  Additionally,  information 
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explaining  the  effect  that  other  seal  parameters  have  on  rotor  stability 
is  needed.  For  example,  increasing  seal  diameter  is  expected  to  enhance 
stability  because  increasing  rotor  and  seal  surface  area  should  increase 
fluid  friction  —  an  energy  dissipating  and,  therefore,  stabilizing 
phenomenon.  Similarly,  increasing  seal  strip  separation  in  Figure 
1.6)  should  also  be  stabilizing. 

If  the  individual  effects  of  such  parameters  can  be  isolated, 
then  design  criteria  will  have  a  basis  for  development.  Such  information 
can  allow  derivation  of  formulas  relating  seal  parameters  and  rotor 
stability. 

The  purpose  of  this  thesis  was  to  provide  such  information. 
Analytical  verification  of  stability  threshold  increase  due  to  divergent 
seal  geometry  was  its  principle  aim.  Since  stability  was  defined  in 

Disc 

Left  Bearing  Right  Bearing 


Figure  1.8  -  Myrick  Rotorbearing  System  Model  with 
Single  Labyrinth  Seal 
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terms  of  frequency  the  analytical  model  developed  xn  this  thesis  produced 
root  locus  output.  This  model  used  a  rotor  bearing  system  similar  to  one 
described  by  Myrick,  adding  a  single  labyrinth  seal  (Figure  1.8).  It  was 
also  similar  to  Wright's  experimental  apparatus.  It  simulated  divergent, 
convergent,  and  straight  seal  geometries  as  well  as  variations  in  seal 
dimensions,  pressure,  and  rotor  speed.  Stabilizing  or  destabilizing 
trends  were  available  after  each  parameter  was  separately  incremented. 

This  thesis  presentation  will  first  describe  its  model 
physically.  Components  will  be  identified  and  located  relative  to  each 
other.  Operating  conditions  such  as  pressure  and  temperature  will  be 
presented . 

Derivation  of  the  analytical  model  will  follow  system 
description.  This  chapter  will  present  assumptions,  derive  equations  of 
motion,  and  prepare  these  equations  for  computer  programming.  The 
schemes  needed  for  such  programming,  as  well  as  a  listing  of  model 
parameter  exercises,  will  also  be  presented.  General  information  needed 
to  understand  these  computer  programs  will  be  made  available. 

Test  results  will  present  and  explain  output  generated  by  these 
computer  programs.  Eigenvalues  will  be  graphed  on  root  locus  plots. 
Trends  will  be  identified  and  explained. 

Finally,  conclusions  based  upon  test  results  will  be  presented. 
This  final  chapter  will,  ultimately,  satisfy  the  aim  of  this  thesis.  It 
will  state  its  success  or  failure. 


Chapter 


SYSTEM  DESCRIPTION 


2. A  Rotor  Bearing  System 

The  rotor  bearing  system  described  by  Myrick  served  as  the 
basis  of  this  study  (Figure  1.8).  It  consisted  of  3  masses  (2  journals 
and  a  disc)  separated  by  2  flexible  shafts.  Myrick’ s  purpose  was  the 
study  of  rotor  bearing  system  whirl  using  realistic  incompressible  film 
hydrodynamic  journal  bearings.  Since  analytical  description  of  such 
bearings  was  complex,  lengthy,  and  required  considerable  computer 
simulation  time,  Myrick 's  bearing  model  was  modified  to  assure  an 
inf initessimally  short  bearing,  i.e.,  a  bearing  with  fluid  flow  in  the 
circumferential  direction  only  (Figure  2.1)  described  by  the  following 
reduced  version  of  the  Reynold's  equation: 


)z  \  i?/  r 

Here,  ©  was  an  independent  variable  representing  circumferential 
angular  distance.  Its  datum,  per  Figure  2.1,  coincided  with  the  line  of 
centers  at  the  position  of  maximum  journal  to  bearing  clearance.  Minimum 
journal  to  bearing  clearance  was  represented  by  h.  Dynamic  viscosity  was 
represented  by/c>  journal  surface  speed  by  U,  and  bearing  fluid  pressure 
by  P.  Axial  distance  was  represented  by  z. 

This  short  bearing  simplification  is  supported  by  Myrick' s 
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conclusion  that  hydrodynamic  moments  in  finite  length  journal  bearings 


\ 


Figure  2.1  -  Short  Bearing  Circumferential  Cross-section 


are  small  compared  with  shaft  bending  moments.  Additional 

approximizations  consistent  with  this  short  bearing  assumption  are:  [“] 

(1)  Bearing  and  journal  curvatures  are  small  compared  with  film 
thickness. 

(2)  The  pressure  and  journal  to  bearing  clearance  (h)  gradients  in  the 
axial  direction  are  much  less  than  those  in  the  vertical  direction  (y  in 
Figure  2.1).  i.e. , 
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If.  <<  A£  -lib-  <<  ik. 

*2  dj  )2  ^ 


(Short  Bearing) 


(3)  h  =  c+e  cos@  =  minimum  journal  to  bearing  clearance  where  c  equals 
radial  clearance  between  a  concentric  journal  and  bearing. 

(4)  The  bearing  fluid  is  incompressible  and  has  constant  viscosity. 

(5)  Bearings  are  rigid  and  stationary. 

The  two  flexible  shafts  were  assumed  to  be  massless.  They  had 
circular  cross  sections,  homogeneous  compositions,  and  no  damping 
properties  (completely  elastic).  These  assumptions  were  consistent  with 
finite  element  modeling  necessary  for  computer  simulation.  Shaft 
stiffness  in  bending  was  described  by  the  relation: 


(2.2) 


The  centrally  located  disc  and  both  journals  were  assumed  to 
behave  as  point  masses.  i.e.,  tilting  and/or  deformation  was  neglected. 
Furthermore,  axial  displacement  was  also  neglected.  Because  eccentricity 
ratios  (£  *  e/c)  in  this  analysis  did  not  exceed  .5,  Myrick's  assumption 
of  constant  rotor  speed  (u>)  applied  for  this  model.*  This  implied  an 
absense  of  shaft  twist  since  constant  rotational  speed  for  all  masses 
required  identical  rotational  speed  for  all  masses.  Therefore,  elastic 
shaft  deformation  occured  in  bending  only  and  mass  displacement  occured 
only  in  a  plane  transverse  to  its  axial  steady  state  rotor  position. 


*  Myrick's  upper  limit  for  constant  rotor  speed  was  £.  ■  .7 


2.B  Labyrinth  Seal  System 


The  model  labyrinth  seal  (Figure  2.2)  was  concentrically 
located  about  the  central  disc  for  any  and  all  rotor  bearing  system 
steady  running  conditions  (constant  speed) .  Air  was  used  as  the  working 
fluid  through  and  within  the  seal  to  facilitate  possible  experimental 
duplication  of  this  analysis.  Downstream  pressure,  P  ,  was  constant  and 

D 

identical  for  most  simulations.  It  represented  exhaust  into  an  ambient 
environment .  In  accordance  with  Kearton  DO  and  Kostyuk  DO  ,  fluid 
temperature  was  considered  constant  and  axial  fluid  kinetic  energy  was 
considered  to  be  completely  destroyed  within  the  seal  chamber. 


Figure  2.2  -  Model  Labyrinth  Seal 


Chapter  3  ANALYTICAL  MODELING 

3. A  Overview 

Analytical  modeling  of  this  system  was  actually  accomplished  in 
three  phases.  Each  phase  consisted  of  equation  derivation  and  computer 
simulation  for  a  particular  analytical  model.  Equation  derivation  for 
all  models  occured  first  followed  by  explanations  of  simulation 
procedures . 

Phase  1  developed  the  nonlinear  rotor  bearing  system  model. 
Equations  of  motion  were  derived  for  the  rotor  bearing  and  labyrinth  seal 
systems  separately.  Using  a  finite  difference  time  integration  scheme, 
rotor  bearing  system  operation  was  compared  with  Myrick’s  observation  of 
stability  threshold. 

In  phase  2,  these  rotor  bearing  system  equations  were 
linearized,  reduced  to  first  order  form,  and  assembled  in  matrix  form  [lj*. 

0  X 

-M'V  -fcf'jB 

This  formulation  allowed  the  eigenvalues,  ^  ,  to  be  determined  using  the 
QR  algorithm.*  Corresponding  eigenvectors  were  then  obtained  by 
subjecting  the  matrix  in  equation  3.2  to  Gaussian  Elimination. 

*  QRHMOD  -  available  on  METAPE  and  DYNSYS  Mechanical  Engineering 
Department  computer  files. 


-  > 


% 


(3.1) 
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(3.2) 


Eigenvalues  allowed  root  locus  system  representation  while  eigenvectors 
described  relative  state  variable  magnitudes  in  each  discrete  mode  of 
system  behavior.  When  results  from  these  formulations  matched  those  of 
Shapiro  and  Colsher  £lj  rotor  bearing  system  modeling  verification  was 
confirmed  since  they  similarly  modeled  Myrick's  system*.  With  this 
assurance,  a  labyrinth  seal  model  was  incorporated  into  this  rotor 
bearing  system  model  during  the  third  phase. 

Labyrinth  seal  nonlinear  modeling  was  more  difficult  than  that 
for  the  rotor  bearing  system  since  both  time  and  circumferential  seal 
distance  were  bases  for  state  variable  integration.  It  required  the  use 
of  a  finite  element  labyrinth  seal  model,  needed  to  effect  seal 
circumferential  distance  integration.  The  expense  of  such  an  effort,  in 
computer  simulation  time,  was  unjustified.  Contrary  to  experience  with 
the  rotor  bearing  system  nonlinear  modeling,  comparison  with  a 
nonexistent  and  proven  seal  model  was  impossible.  Comparison  of 
nonlinear  and  linear  seal  models'  results  had  little  value  since  both 
models  derive  from  the  same  equations.  Therefore,  nonlinear  labyrinth 
seal  modeling  did  not  occur  beyond  derivation  of  system  equations  as 
required  for  linear  modeling. 

Linear  labyrinth  seal  modeling  took  a  form  somewhat  different 

*  Calculation  details  are  not  available.  Figures  4.4b  and  4.5b  contain 


these  results. 
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from  that  for  the  rotor  bearing  system.  Since  linearization  of  equations 
of  motion  could  not  eliminate  partial  derivative  terms,  assumed  solution 
forms  were  substituted  for  each  state  variable.  The  result  was  four 
first  order,  time  dependent  differential  equations  with  solution  form 
coefficients  as  dependent  variables.  The  resulting  matrix  was  subjected 
to  the  QR  algorithm  and  Gaussian  Elimination  both  separately  and  coupled 
to  the  rotor  bearing  system's  linear  model  matrix.  As  a  result, 
eigenvalues  and  eigenvectors  from  coupled  and  isolated  systems  were 
compared  for  a  complete  stability  analysis. 


3.B  Rotor  Bearing  System  Modeling 

3.B.1  Nonlinear  Equation  Derivation 

Initial  analytical  modeling  of  Myrick's  rotor  bearing  system 
required  adequate  description  of  bearing  film  forces  on  each  journal. 
Beaman  jjLjfJ  has  accomplished  such  an  analysis. 

Starting  from  the  reduced  Reynold's  equation  (equation  2.1),  an 
expression  describing  bearing  pressure  per  unit  axial  distance  was 
obtained: 


F  ~  ~  AA^ 

5L 


U  Ik 

rW  >9 


4- 


ail 

h’ 


(3.3) 


where  r  ■  journal  radius 

Forces  normal  to  and  in  the  eccentricity  direction  (e  -  Figure  3.1)  were 
derived  by  integrating  bearing  film  pressure  around  the  journal 
circumference: 


_  .  'Scr  .i 


22 


These  forces  are  generated  by  the  bearing  film  and  act  on  ea<h 
journal's  centroid.  After  integration,  equations  3. A  become: 


(i -cY  (1-tT* 


(3.3a) 


Xc 


hlTg.  _  ?£'g. 

a(i-€.*yA‘  (/-s*)1 


(3.5b) 


/  ■  I  1  y-  1 

(note:  £  and  <f>  involve  differentiation  with  respect  to  ) 

Conversion  to  a  standard  x-y  coordinate  system  (Figure  3.1)  is  easily 
accomplished  using  the  following  relations: 


Fx  =  -S  ir\  <f>  +  (F*)cos<|> 


(3.6a) 


Fy  “  (  F^  cos  4-  (F^)  s  m  <t> 


(3.6b) 
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Equations  3.6  apply  to  both  left  and  right  bearings.  They  can  be 
considered  to  be  generalized  forces  acting  on  each  journal's  _entroid 
(Figure  3.2)  £lxQ  • 

Figure  3.2  depicts  Myrick's  model  with  nonlinear  short  bearing 
forces  and  linear  x  and  y  direction  springs  connecting  each  journal  to 
the  central  disc.  From  this  figure  and  Chapter  2  assumptions,  two 
generalized  coordinates  for  each  mass  were  assigned  (Figure  3.3). 

Since  all  constraints  imposed  upon  this  system  were  holonomic, 
Lagrange's  equations  were  used  to  derive  its  equations  of  motion  [**]• 
Starting  with  the  central  disc,  a  kinetic  coenergy  expression  was 
derived. 

T4*-  -ir  fY)d  (x,z4-  X?)  4-i  (J~ 

Disc  potential  energy  was  stated  as: 

Vj  -  -t  k*  Cx,  ~  ^5^  "J 

-f  “kk>j  +  (x9.~xJn) 

Combining  these  expressions,  the  disc  Lagrangian  was  expressed: 

1,=  T/-Vj  =  + 

— kkx'|x.-x,')*+Cx,-x^]--tkr 

Lagrange's  equations  (equation  3.10)  for  the  disc  were  then  formulated, 
realizing  that  all  forces  acting  on  the  disc  were  conservative  (S’  -  0) 
and  that  rotor  angular  velocity, W  ,  was  constant  and,  therefore,  not  a 


(3.7) 


(3.8) 


generalized  coordinate. 


d_  fix\  _  a- 
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(3.10) 


For  coordinate  X^  (Figure  3.3  -  G.C.  assignment): 

X|  -  —  k>(  *s) 

ma 

Similarly,  for  coordinate  X2: 


X*  --  K,  {Xl  -  X,)  _ 

ma 


(3.11) 


(3.12) 


Equations  of  motion  for  each  journal  were  similarly 
because  bearing  film  forces  include  nonconservative 
represented  as  Generalized  Forces,  i.e.. 


derived.  However, 
damping  they  were 

(3.13a) 


f V 


(3.13b) 


Proceeding  for  the  left  journal  as  in  the  derivation  of  disc  equations: 

Tr*=-k«rU>*;)+iIMwl  o.  1 


V]~L  “  ( X3  A  ^  ^  (3. 15) 

^TL  “  +  *+  )"t"  ^-AT  ^  C-Vs 

~  ^  (*•<  ”  *x)~\ 


t 


(3.16) 


Employing  Lagrange's  Equations  (equation  3.10)  and 
equations  3.13. 


*,  =  -Jfj.  -*) 


YAT 


-Xi)  4. 

rvoT  ^  r 


Similarly,  for  the  right  journal: 

Vr*  X,Y  +  i  kj  Of i-X’aY' 


X5  —  for  lVs  ~  4-  F> cr 

i^x  ™  x 


incorporating 


(3.17) 


(3.18) 


(3.19) 


(3.20) 


(3.21) 


^  L  =  foj  iX|,_  4 


(3.22) 
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Equations  3.11,  12,  17,  18,  21,  and  22  were  integrated  with 
respect  to  time  using  a  Runge  Kutta  finite  difference  technique.  To  do 
this,  it  was  convenient  to  reduce  these  2nd  order  differential  equations 
to  1st  order  form.  Defining  the  following  additional  variables  allowed 
such  reduction: 


*7  -  (3.23a) 

X?  ~  X a  (3.23b) 

-  Xj  (3.23c) 

> 

Xxb  '  XH  (3.23d) 

^  Xs  (3-23e> 

X,x  -  XL  (3. 23f) 


Since  these  variables  were  not  independent  of  x,  through  x,  but  still 

l  o 

specified  system  behavior,  they  were  labeled  as  state  variables  [,<]. 
Rotor  bearing  system  equations  of  motion  were  then  restated: 


X-j  ~  k *  " X$) 

W4  "U 


(3.24a) 


(3.24b) 
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(3.24c) 
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3.B.2  Linearisation 


Since  equations  3.24c  through  f  have  nonlinear  components  they 
were  not  suited  for  a  frequency  domain  solution,  as  required  for  root 
locus  output.  Therefore,  these  bearing  film  force  terms  had  to  be 
linearized. 

Beaman  had  linearized  equations  3.5  to  achieve  the 
following  form: 


3* 
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(3.25) 


where ie'  and  Sf’  indicated  differentiation  with  respect  to  TT *  tu/r 
and:* 


*  Stiffness  and  damping  coefficients  appear  in  Appendix  B. 
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(3.26a) 


(3.26b) 


Equation  3.25  was  converted  to  the  Figure  3.1  coordinate  system  using  the 
following  transformation: 


X  -  e.  sm<t> 


(3.27a) 


'j  =  -  e_  cos  <(> 


(3.27b) 


The  following  transformed  solution  resulted: 


where,  after  equations  3.26  were  incorporated:* 


-u 

c 


(3.29a) 


*  Stiffness  and  damping  coefficients  appear  in 


Appendix  B. 
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(3.29b) 


Bearing  force  terms  in  equations  3.24c  through  f  were  then  expressed  in 
linear  form  using  the  equation  3.28  revised  bearing  coefficients. 


— 
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(3.30a) 


(3.30b) 


(3.30c) 


(3. 30d) 


With  these  substitutions,  equations  3.23  and  3.24  became  linear  rotor 
bearing  system  equations  of  motion.  They  were  then  arranged  per 
equations  3 . 1  and  3.2. 

The  square  matrix  in  equation  3.1  is  an  algebraic  manipulation 


of  the  more  fundamental  linear  homogeneous  matrix  equation  of  motion 
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(3.31a) 


Where,  a  represents  a  generalized  coordinate  column  vector  and  in 

equation  3.1,  represents  a  state  variable  column  vector  (not  all 
variables  are  independent) : 
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(3.31b) 


Equation  3.3la  is  reduced  to  1st  order  form  by  doubling  the  number  of 
variables  as  in  equation  3.23.  Equation  3.1  is  further  rearranged  to 
facilitate  solution: 


-  A](x|  --  o 


where: 

A  = 


o 


and,  [l9] 
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(3.32a) 


(3.32b) 


(3.32c) 


Equation  3.32a  was  solved  for  ,  the  system  eigenvalues. 
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These  eigenvalues  were  then  used  in  the  following  formulation  to 
determine  the  system  eigenvectors  (  v  )  corresponding  to  each  eigenvalue: 


(3.33) 


A  was  then  developed  from  equations  3.23,  3.24,  and  3.32b: 
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(assume  blank  spaces  contain  zeros) 


(3.34) 


3.C  Labyrinth  Seal  Modeling 


3.C.1  Derivation  of  Equations  of  Motion 

Figure  3.4  depicts  the  model  labyrinth  seal  surrounding  the 
central  rotor  bearing  system  disc.  In  this  analysis  the  following 
additional  assumptions  were  required: 


(Note:  h  now  represents  minimum  seal  base  to  disc  radial  clearance) 
Figure  3.4  -  Single  Labyrinth  Seal  Model 

1)  Working  fluid  (air)  behaves  as  an  ideal  gas:  H 
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(3.35) 

P  =  ^  RT 

(^  =  Fluid  Density,  R  =  Gas  Constant) 

2)  The  disc  is  concentric  within  the  rigid  and  stationary  seal  during 
steady  running. 

3)  Fluid  axial  kinetic  energy  is  completely  destroyed  within  the  seal 
chamber.  5] 

4)  This  flow  process  is  adiabatic 

5)  Fluid  temperature  (T)  is  constant  throughout  the  seal  KJ. 

6)  Circumferential  fluid  velocity  during  steady  running  is  h  that  of 
disc  surface  speed  and  is  in  the  direction  of  disc  rotation:  (Figure 
3.5). 

-  rvO  (3.36) 

A 


3 


where, 

K,  =  -otnt 

n 

K,  -  .0m3 

n  •l.H’ 


1 0  Re.  ^  I  A 

10s  <  Kt  <  \0' 


8)  Flows  through  seal  strips  A  and  B  (Figure  3.4)  are  each  modeled  as 


that  of  an  isentropic  nozzle  affected  by  a  discharge  coefficient jj.2] 

9)  Fluid  dynamic  viscosity  is  assumed  constant. 

10)  Axial  flow  is  always  positive,  i.e.,  PA^  Pfi  (constant  pressures), 


®  ©  © 


r 


s 


Section  A-A 

“  seal  angular  distance  (independent 

variable) 

*  seal  circumferential  distance^. 


Figure  3.5  -  Labyrinth  Seal  Control  Volume 
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As  in  any  compressible  fluid  mechanics  analysis,  the  concepts 
of  continuity,  momentum,  state,  and  previous  assumptions  were  applied  to 
the  control  volume  of  Figure  3.5.  Because  of  the  adiabatic  flow 
assumption  an  energy  balance  was  not  necessary.  One  differential 
equation  for  each  generalized  coordinate  resulted.  In  this  case,  two 
equations  were  expected  corresponding  to  the  generalized  coordinates  of 
chamber  pressure  (P)  and  circumferential  fluid  velocity  (u) .  Minimum 
disc  to  seal  base  clearance  (h) ,  although  variable,  was  not  independent 
since  it  was  directly  related  to  disc  behavior. 

Equation  3.38  expressed  the  mass  flow  balance  or  continuity  of 
the  system: 

-LU^  4  i  (pMm)  -  -  'V'j  (3.38) 

It  is  ~  [l-Cont] 

Considering  shear  forces  on  circumferential  and  radial  interior  (within 
chamber)  seal  surfaces  only,  a  momentum  balance  was  performed: 

(?h  ~  -  r  jJs  +  tJi+akV* 

(3.39a) 

£  2-Mon^ 

it 

where  friction  stresses  were  represented  by: 

I(  -  i  ft  ?  (rw-u) 


(3.39b) 
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(3.39c) 


Differentiating  equation  3.38  and  substituting  equation  3.35 


yielded: 


Ik  vr  r  (».-*,')  -  IP  3k. 

RT  it  RT  Jt 


f/Rk 

-1-  U^k\ 

1 - 

M 

3 

+ 

%T 

l  is 

^  J 

(3. AO) 


Equation  3.39a  was  also  differentiated.  After  substituting 
equations  3.39b  and  c  it  became  equation  3.  A  It, 


Ifk  3m  -  _J> 
RT  it 


? 3R  k  3  P 

i  5 


■f.  PuVC* 

3a  KT  it 


+  uh  _^P  +  UZf] 

)s 


a-RT 

Tf  u  2k  +  +  a  u K  ^ 

V  <)t  is  is 


(3.A1) 

|~2-MonTj 


Equation  3.  AO  required  the  following  mass  flow  rate 
substitutions 


%=  0*-^anr/k7  & 


RT 


Ma 

-I \k.*A  ("*•>/»<*- 'll 
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43.42) 


( i * (^)hn 


m6-(K~k);LTrr/k~"  PB 


(3. A3) 


Where  CWA  and  represented  empirical  discharge  coefficients  for  seal 
strips  A  and  B  respectively: 


CyVA  —  f  (  v/p^)  4-  Constant 


(3.44a) 


Cw0  -  -f  (? /ft)  +  constant 


(3.44b) 


Mach  numbers  and  were  functions  of  chamber  pressure,  P.  (see 
equations  #.  *3  9>2f) 

3.C.2  Linearization 


t 


r-  p.  +  sp 


K  '  hn  + 


U  - 


f i  "  fin  +• 

Jo  “  fon  +  f  fo 


™A«  + 

m6=  w4n  4-Sf^ 

MA  ~  /^An  +  SM* 
—  M0„  +  SM<3 


C-WAm  ■*“  ^  ^WA 


+•  J  ^W3 


=  Chamber  pressure 


=  Seal  Base  to  Disc  Clearance 


=  Circumferential  fluid  velocity 


=  Friction  Factor  at  Disc  surface 


=  Friction  Factor  at  Seal  Base 


=  Mass  Flow  Rate  at  Seal  Strip  A 


=  Mass  Flow  Rate  at  Seal  Strip  B 


=  Mach  number  at  Seal  Strip  A 


*  Mach  number  at  Seal  Strip  B 


•  Discharge  Coefficient  at  A 


■  Discharge  Coefficient  at  B 
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(3.45a) 

(3.45b) 

(3.45c) 

(3.45d) 

(3.45e) 

(3 . 45f ) 

(3.45g) 

(3.45h) 

(3.451) 

(3.45J) 

(3.45k) 


Mass  flow  rates  in  equation  3.40  were  linearized.  Therefore,  MA  and  Mg 
required  definition.  Since  these  were  axial  flow  rates  only,  they  were 
not  functions  of  circumferential  velocity,  u.  The  following  general 
relations  applied  and  were  substituted  into  equations  3.45f  and  g: 
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SP  + 


SP  + 


St* 


The  "K"  coefficients  in  equations  3.46  were  derived  from 
equations  3.42  and  3.43  with  the  following  substitutions 


K 


0  + 


M0v/li-  (*¥)MeX 


K,  --  /L  Pa 

v  pr 

K-  [J_  PB 

V  RT 


Cy,A  =  +  -1/72372? 

Cv,B  =  -.Jliwm  +  .7172  3  72? 


(3.46a) 


(3.46b) 


linearization  of 


(3.47a) 


(3.47b) 


(3.47c) 


(3.47d) 


(3. 47e) 


(3 . 47f ) 


A.«  =  U  -»)  ifrr 


(3.47g) 


A„„  --  (K-t'j  a-rrr  ,  «•«»> 

J 

Equations  3.47e  and  f  represented  linearized  approximations  of  J.A. 
Perry's  graphical  discharge  coefficient  relation  for  sharp  edged  orifice 
meters.*  Incorporating  equations  3.47  into  equations  3.46  yielded 
equations  3.48. 

-  K  {  +  * K\  S  Cw<]  (3  48a) 

+  r/u] 

fm„  =  K'b^A,„Jcw<„  +  T(Me\  ^C«.0]  o.48b) 

Previously  undefined  terms  in  equations  3.48  required  definition:** 


-  3-TrrSlA  j  ABri  =  (k-^atr r 


*  This  graph  appears  on  page  100  of 

**  "K"  coefficient  expressions  appear  in  Appendix  B. 


=  ^irrk  j  Ab„-  (k„-  t>)  Stt  r 
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(3.49b) 


Cwa^  -  -.virtu  ^  + 

(3.49c) 

c*.,  -  .inawi 

(3.49d) 

£Cvv*  "  ^gCWA  ^  ^ 

(3.49e) 

5  Cv*fi  ~  ^jcwj  ^  ^ 

(3. 49f) 

'  K^Ma 

(3.49g) 

J(5*C^o5)  r  ICsj-mo  ^  ^3 

(3.49h) 

S  Ma  -  ^  ^ 

(3.491) 

*Jmb  SP 


(3.49j) 


-  f  OC) 


=*  Cm^)„  - 


Therefore; 

£  =  PjfA fA  “  PjfAIAT  ^  P 

S  (?(/*V^\  -  Pj^mo  ai 5  ^STtfbr 


Coefficients  in  equations  3.46  were  then  described  in 


equations 

3.47,  3 

.48,  and  3.491 

Pp*  A  ~ 

Pa 

MW 

^An  (CwAr\  Pff/HAr 

+  #*A  *^1 

o?77T 

L 

Pp*»B 

P* 

+■  ^  Wn  Pjcwe^l 

a*rrr 

Kkm*  =  CwAn  { 


(3.49k) 


(3.491) 


(3.49m) 


(3.49n) 


terms  of 


(3.50a) 


(3.50b) 


(3.50c) 
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(3. 50d) 


Realizing  that,  during  steady  motion,  nominal  axial  flow  rates  at  A  and  B 

(Figure  3.5)  were  equal  (mt  =  m_  ) ,  equations  3.40  and  3.41  were 

An  an 

expressed  in  terms  of  the  small  motion  substitutions  of  equations  3.45: 
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The  friction  related  coefficients  in  equation  3.52  exist  because  friction 


factors  f.  f  are  functions  of  state  variables  -  p,  u,  and  h.* 

1,0  r 


(3.53) 


Equations  3.51  and  3.52  were  then  expressed  more  compactly! 

~  4.  ^4-  [Kik]^  4-  jCtjUk 

3  "t  3  tr  ^ 


Mai*  + 1 
L  1 

^  1  p j|  a  s  p  ~  ^ 

“  J  ^  5 

Q-Cont] 


*  Appendix  B 
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These  linearized  continuity  and  momentum  equations  were  not 
compatible  with  the  linearized  rotor  bearing  system  model  (equation 
3.32c)  because  partial  derivations  with  respect  to  both  time  and 
circumferential  distance  existed.  The  linear  rotor  bearing  system  model 
was  based  upon  time  as  the  only  independent  variable.  Consequently, 
these  linearized  seal  equations  could  not  be  coupled  within  the  A  matrix 
of  equation  3.34.  Therefore,  it  was  necessary  to  reconfigure  equations 
3.54  and  3.55  such  that  partial  derivatives  no  longer  appeared  and  time 
became  the  only  independent  variable.  One  approach  involved  substitution 
of  the  following  assumed  solution  forms  into  equations  3.54  and  3.551 

$P  -  (Ap  CosV  4-  Bp sin  f)  o.56a) 

§u  -  +  Bi*  Sir  Vs)  (3.56b) 

Sk  -  tU  ^  (3.56c) 


Refering  to  Figure  3.4  (Seal  Model)  and  remembering  that,  during  steady 


running. 

the 

disc 

is  concentric  within  the 

seal,  equation  3.56c  was 

rewritten 

: 

J  =  0  a*<$ 

s  k  =  - 

cesV 

<f>*  -  0 

II 

WO 

cos^l 

V* 

(3.57) 

Sk  = 

-H 

toil 

rckt 

«i 


Therefore,  H  represents  disc  displacement  in  the  y  direction  (Figure 
3.4),  i. e. , 


^  =  Yt* 


(3.58a) 


and,  from  equation  3.56c: 


=  "  H  COS  If 


(3.58b) 


therefore, 


Sk 


(3.58c) 


or 


> 


h  =  y  a.a  h  =  y 


(3. 58d) 


Substituting  equations  3.56  into  equations  3.54  and  3.55  yields  time 
dependent  differential  equations  without  partial  derivatives; 
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f~  ^1-P^  +hjkaL  NAu+Z ^IkAH  (3‘59b) 

'  r  k’lPt  /  \  ^irt  J  V  PK'ipt-  /  vKipJ 


Am  ~/*"K»p  —  K&pt  ^ip  \  Ap  -j-  /~ K^pt-  K^pj  —  K~af5 

Vr  K-tPtr  *r»ut.  ^  *; 


B, 


K^Ut-  ^lPf  ^4«t 


XMt 


-f/’  Au  4-  (~  ^zpt-  Kjus  —  ftruj 

'  KsLut  /  '  P 


B 


^tPt-  ^iut 


r  K’vut- 


(3.59c) 


YJ^ifL^tL  ■+  ^  +  Z  fcfr  ^at  4-  fakt- 

V  ^ire  ^*.ut  K’zuwt]  »  ^lffr  K***t 


w 


in 

^iut 


*4-  + 

\r  Kipr 


4  /~^tP  —  ^iPt 
V  K»wt  ^ipc 

-fcwA  Bu 

^Jut  1 


4" 


fcft  - 

r  tCifr  fcrVut- 


(3.5vi, 


Linearized  seal  equations  were  then  present  in  a  form  similar 


to  equation  3.32c: 


Because  H«=Y  represented  disc  displacement  in  the  y  direction  it  also 

corresponded  to  X„.  Consequently,  H  *  XQ.  Equation  3.60,  therefore, 
Z  o 

included  the  effect  that  the  disc  has  upon  the  seal.  To  complete 
coupling  between  seal  and  rotor  bearing  systems,  seal  force  exerted  on 
the  disc  was  also  specified. 

Since  linearization  involved  inf initessimal  variations  of  state 
variables,  seal  force  acting  on  the  disc  was  so  represented.  (See  Figure 
3.4  -  Seal  Model) 
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r3TT  f  a>sn  Hr  a  r 


(3.61b) 


To  accommodate  these  integrations,  equation  3.56a  can  be  expressed  as  in 
equation  3.62! 


f  -  Ap  (i)  y  +  Bp  (t)  SI»  Y 


(3.62) 


Substitution  into  equations  3.61  yields: 


fF„  =  Bp  -rr  rJl 


(3.63a) 


JFj  -  Apirri 


(3.63b) 


Since  the  datum  for  independent  variable  time  in  equation  3.32c  was 


arbitrary,  t 

Ar  (*.')  = 


Br  ^  = 


=  0  was  chosen. 


Consequently, 


and : 


(3.64a) 


(3.64b) 


Therefore,  seal  to  rotor  bearing  system  coupling  was  accomplished. 
Equation  3.65  depicted  the  A  matrix  of  equation  3.34  expanded  to  16  x  16 
to  include  equations  3.60  and  3.63! 


(Blank  spaces  contain  zeros:  upper  left  hand  corner  is  the  A  matrix  of 
equation  3.34) 

(3.65) 

Equation  3.65  was  then  used  to  generate  eigenvalues  and 
eigenvectors  for  the  rotor  bearing  and  labyrinth  seal  systems,  both 
isolated  and  coupled. 

3.D  Simulation  Procedure 

This  section  will  explain  the  programming  schemes*  needed  to 


*  Listed  in  Appendix  A. 
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exercise  the  nonlinear  rotor  bearing  system  (equations  3.23  and  3.24)  and 
the  A  matricies  of  equations  3.34  and  3.65.  Programming  of  the  nonlinear 
rotor  bearing  system  model  will  be  discussed  first.  Programming  of  both 
rotor  bearing  and  labyrinth  seal  systems  linear  models  will  be  explained 
next  to  be  followed  by  a  description  of  parameter  variations. 

The  parameter  variations  section  of  this  chapter  will  focus 
primarily  upon  labyrinth  seal  physical  variables.  Myrick’s  rotor  bearing 
system  will  experience  minor  changes  only,  to  minimize  impact  on  possible 
follow  on  experimental  verification.  Such  testing  will  provide  data  for 
analysis  in  Chapter  4. 

3.D.  1  Isolated  Rotor  Bearing  System  Testing 

3.D.1.A  Nonlinear  Model  Programming  (RW5) 

As  previously  mentioned,  nonlinear  rotor  bearing  system 
modeling  was  intended  as  a  modeling  accuracy  check.  Its  basis  was  a  4th 
order  Runge  Kutta  finite  difference  integration  of  equations  3.23  and 
3.24.  Although  €.  and  ^  in  Figure  3.1  (Bearing  description)  were 
directly  related  to  bearing  generalized  coordinates  X_  through  X,,  they 

J  D 

were  also  separately  integrated  to  simplify  bearing  force  calculations 
(equations  3.5).  The  coordinate  assignments  of  equations  3.23  and  3.24 
did  not  change.  Additional  bearing  state  variable  assignments  were  per 
equations  3.66.* 

X\y  S  4*  '')  *  Left  Journal  (3.66a) 


Eccentricity  ratio 


r 


where  c  =  Journal  to  bearing  radial  clearance 


bk 


X,H 


=  Left  Journal 

Attitude  Angle  (3.66b) 


=  Right  Journal  (3.66c) 

Eccentricity  ratio 


&r) 


=  Right  Journal 

Attitude  Angle  (3.66d) 


First  derivative  expressions  for  these  new  state  variables  were  also 
required: 


(3.66e) 


(3. 66f ) 
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( X*  K  ±  -X,  O  =  i 


^15  C 


(3. 66g) 


F,t  =  -  ( x,  x„  -  O  =  4>r  <3-66h> 

Since  journal  surface  speed  (U  =  ru>)  in  equations  3.5  required 
specif ication  of  angular  velocity  an  expression  from  Cameron  0°]  was 
used : 


fWc*(  1- 

a 

-7T>c  Jpcr  f 

V  + 1 ! 

1  'lx 

1 

1 

(3.67) 

(W  =  \  total  rotor  bearing  system  weight) 

Journal  length,  r  =  Journal  Radius) 

Equation  3.67  was  compatible  with  bearing  assumptions  made  in  Chapter  2 
and  applied  to  all  rotor  bearing  system  mass  elements.*  Cameron  had  also 
derived  an  expression  for  journal  attitude  angle,  ^  : 


*  Constant  speed  assumed  in  Chapter  2. 
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Since  both  left  and  right  journals  (Figure  3.2  -  Rotor  Bearing  System 
Model)  were  physically  identical  and  since  their  angular  velocities  were 
also  equal,  £•  and  $  were  expected  to  be  identical  for  both  journals, 
i.e.,  symmetry  should  prescribe  identical  behavior  for  both  journals. 
For  verification,  however,  angular  velocity  was  defined  for  each  journal. 

As  in  any  time  dependent  simulation,  initial  conditions  were 
prescribed.  These  were  defined  in  terms  of  an  assumed  journal 
eccentricity  ratio,  £.  ,  and  attitude  angle,  ^  —  the  same  values  at  both 
bearings.  Disc  initial  circumferential  position  was  assumed  equal  to 
that  of  both  journals.  Since  equations  3.67  and  3.68  were  defined  for 
steady  running  only*,  the  initial  value  for  ^  did  not  satisy  equation 
3.68.  Otherwise,  mass  displacement  and,  therefore,  potential  system 
instability  could  not  be  observed.  Instead,  steady  running  conditions 
were  defined  separately  and  differently  with  respect  to  these  initial 
conditions . 

Steady  running  conditions  were  defined  in  terms  of  an  assumed 
eccentricity  ratio  and  corresponding  journal  angular  velocity  (equation 
3.67)  and  attitude  angle  (equation  3.68).  Steady  journal  positions, 
therefore,  were  based  upon:  (1)  static  deflection  due  to  rotor  weight 
and  (2)  angular  displacement  ($)  due  to  journal  interaction  with  bearing 
fluid.  Steady  running  disc  position  was  assumed  equal  to  that  of  both 
journals.  Steady  bearing  film  forces  added  in  the  line  of  centers 
direction  (Figure  3.1).  The  vector  sum  of  these  forces  equaled  rotor 
bearing  system  weight.  Steady  running  conditions,  in  effect,  established 
a  temporary  coordinate  reference  frame  (Figure  3.6).  With 


*  Rotor  operation  without  displacement  of  mass  element  centroids. 


steady  conditions  defined,  computer  simulation  could  begin  by  applying 
different  initial  conditions. 


\ 


»  » 

Figure  3.6  -  Fixed  Inertial  and  Steady  Running  (x  ,y  ) 

Coordinate  Reference  Frames 

Simulation  commenced  with  4th  order  Runge  Kutta  finite 
difference  integration  of  equations  3.23  and  3.24.  State  variable 
increments  were  calculated  for  an  assumed  time  step  (At)  and  added  to 
previous  values  to  yield  state  variable  values  at  time,  t  -  tQ  +  n&t. 
Bearing  force  equations  3.5  appeared  twice,  embodied  with  equations  3.66 
—  once  for  each  journal.  These  equations  were  first  calculated  using 
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£  and  ♦  state  variables  converted  to  a  T  =  tu>  integration  basis: 


(3.69a) 


(3.69b) 


g/  and  <f>'  were  later  multiplied  by  tO  to  allow  finite  difference  time 
integration  —  a  Runge  Kutta  requirement!  Bearing  force  equations  3.6 
appeared  next,  again,  once  for  each  journal.  Steady  bearing  film  forces 
were  subtracted  from  each  bearing  force  (equations  3.6)  —  necessary  when 
force  increments  due  to  transient  loading  were  being  calculated.  Rotor 
bearing  system  state  variable  equations  3.23  and  3.24,  therefore, 
contained  state  variable  and  bearing  force  values  relative  to  a 
particular  steady  running  position  (Figure  3.6). 

Model  accuracy  tests  were  performed  by  setting  steady  running 
conditions  that  were  just  above  and  below  (2  simulations)  the  stability 
threshold  predicted  by  Shapiro  and  Colsher  M  *  This  was  accomplished 
by  prescribing  steady  running  eccentricity  ratios  such  thatu)  in  equation 
3.67  was  less  than  93  rev/sec  for  one  simulation  and  greater  than  93 
rev/sec  for  the  other.  Different  initial  conditions  were  applied  and 
time  dependent  state  variable  values  were  monitored  using  a  simple  plot 
routine.* 


*  MINIPLT — See  Mechanical  Engineering  Dept.  Computer  Files. 


■ttmtmmaHiitSaM 
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To  be  compatible  with  Shapiro  and  Colsher,  testing  below  the 
stability  threshold  should  yield  slightly  damped  state  variable  response. 
Testing  above  the  stability  threshold  should  yield  ever  increasing 
(exploding)  state  variable  values  with  time. 

3.D.1.B  Linear  Model  Programming  (RWE4  ) 

As  mentioned  earlier,  linear  modeling  was  intended  to  produce 
eigenvalues  and  eigenvectors.  Subjecting  the  A  matrix  of  equation  3. 34 
to  the  QRHMOD  routine  accomplished  this  task.  QRHMOD  was  available  by 
calling  METAPE  or  DYNSYS  which  appeared  on  Mechanical  Engineering 
Department  computer  file  tapes. 

To  prepare  equation  3.34  for  QRHMOD,  several  preliminary 
calculations  were  performed.  Beaman's  linearized  bearing 
coefficients  (equations  3.25}  were  calculated  and  converted  to  x-y 
reference  frame  form  (equations  3.21).  Since  linear  analysis  deals  only 
with  infinitessimal  state  variable  variations  from  a  steady  running 
position,  equations  3.67  and  3.68  were  used  to  describe  rotor  speed  (id) 
and  attitude  angle  ( )  respectively,  for  an  assumed  eccentricity  ratio 
(£.).  As  before,  £  values  for  both  left  and  right  journals  were  assumed 
identical.  This  routine  was  used  for  further  comparison  with  Shapiro  and 
Colsher 's  results  and  to  obtain  performance  data  for  the  isolated  rotor 
bearing  system. 

Since  labyrinth  seal  effects  were  the  focus  of  this  thesis  and 
since  changing  the  physical  parameters  of  Myrick's  rotor  bearing  system 
nullified  comparison  of  results,  performance  data,  here,  were  restricted 


60 


to  rotor  speed  variation.  This  was  best  accomplished  by  iterating 
eccentricity  ratio  from  about  £  =  .5  (**>=  45.8  rad/sec)  to  about  £  =  .005 
(°^  =  8745  rad/sec). 

3.D.2  Labyrinth  Seal  Testing 

This  section  discusses  linear  labyrinth  seal  model  programming 
separately  and  in  the  context  of  coupling  to  the  rotor  bearing  system 
linear  model  program  (RWE4) .  Although  explanation  of  isolated  seal  and 
disc  programming  appears  to  be  a  logical  next  step  in  this  discussion, 
this  programming  is  actually  more  conveniently  accomplished  by  modifying 
the  coupled  system  program.  Therefore,  coupled  system  programming  will 
be  discussed  first. 

3-D. 2. A.  Linear  Model  Programming  -  Coupled  Labyrinth 
Seal  and  Rotor  Bearing  System  (RWE4S) 

The  ultimate  intent  of  linear  model  programming  was  to  subject 
the  A  matrix  of  equation  3.65  to  the  QRHMOD  routine.  Since  this  A  matrix 
contained  the  unaltered  rotor  bearing  system  A  matrix  of  equation  3.34  in 
its  upper  left  corner,  expression  of  seal  equation  3.60  coefficients  and 
their  proper  placement  within  equation  3.65  were  accomplished. 

To  do  this,  several  intermediate  tasks  were  performed. 

Initially,  nominal  (steady)  chamber  pressure  (P  )  was  established.  This 

n 

was  accomplished  by  first  assuming  a  Mach  number  at  seal  strip  A  (MAn) 
and  then  solving  for  P^  using  equation  3.70.  This  pressure  (equation 


3.70)  was  then  used  to  calculate  the  nominal  discharge  coefficient  at 
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(3.70) 


seal  strip  A  (C„.  -  equation  3.49c)  and  chamber  pressure.  A  Mach  number 

WAn 

at  seal  strip  B  (M^  )  was  calculated  using  equation  Chamber 

pressure  was  recalculated  using  equation  3.71.  P^  7^  was  compared  with 


PnVil  :=  Aaw  M 

A«r*  C 


■  w  Bn 


k±M flir 

'’■•((i.Nm'V 


(3.71) 


P  „  If  they  differed  by  more  than  .0005xP  _  then  M.  was 

n3.70  n3.70  An 

incremented  accordingly  and  the  process  was  repeated  until  both  P^  values 
satisfied  equation  3.72.*  Since  seal  strips  had  been  assumed  to  act  only 


L,.  I  I -.COOS')  £  Pn,.1|  -  ODOS)  (3.72) 

as  converging  nozzles  both  Mach  numbers  could  not  exceed  1.0.  Mach 
numbers  were  also  made  greater  than  zero  to  avoid  trivial  solutions, 
i. e. : 


0  *  M  -  1 


(3.73) 


Once  P  was  established,  the  axial  flow  coefficients  in 
n 


*.0005  factor  chosen  to  provide  3  digit  accuracy 


equations  3.50  could  be  calculated.  Nominal  circumferential  seal 

fluid  velocity  (un)  was  determined  from  equation  3.36.  Since  steady 

running  disc  position  had  been  assumed  to  be  concentric  within  the  seal, 

nominal  disc  to  seal  base  clearance  (h  )  was  constant  for  all 

n 

simulations.  All  that  remained  was  to  describe  the  coefficients  in 

equation  3.60  and  arrange  them  together  with  rotor  bearing  system 
coefficients  per  equation  3.65. 

3.D.2.B.  Linear  Model  Programming-Isolated  Seal  and  Disc 
(RWE4P) 

With  the  coupled  system  A  matrix  of  equation  3.65  completely 
programmed,  modifications  for  isolated  seal  and  disc  tests  were  then 
described.  All  that  was  necessary  was  to  replace  the  coupled  system  A 
matrix  with  the  seal  PC  matrix,  consisting  of  the  4x4  lower  right  corner 
of  equation  3.65.  In  order  to  properly  reset  this  matrix  between 

simulations,  it  was  filled  with  zeros.  The  calculation  routines  within 
QRHMOD  made  this  necessary  whenever  parameter  variations  required 

repeated  use  of  QRHMOD. 

3.D.3.  Parameter  Variations 

In  order  to  fully  understand  labyrinth  seal  and  rotor  bearing 
system  interaction,  certain  seal  and  disc  parameters  were  varied 
separately.  Corresponding  linear  model  programs,  both  isolated  and 
coupled  (RWE4P  and  RWE4S),  were  used.  Generally,  straight  seal  geometry 


Disc  Speed  (u^)  variation  was  accomplished  by  iterating 
eccentricity  ratio  (£.  )  as  in  RWE4.  Variation  of  seal  chamber  axial 
length  (.?)  was  more  restrictive,  however,  since  too  narrow  a  seal  chamber 


®  ©  ® 


Section .A  -  A 

Figure  3.7  -  Isolated  Labyrinth  Seal  Surrounding  Rotor 
Bearing  System  Disc 

would  violate  the  Chapter  2  assumption  requiring  complete  axial  kinetic 
energy  destruction  within  the  seal  chamber.  To  determine  this  transition 
point,  JL  was  incremented  upward  starting  atJ?“.001  in.  It  did  not  exceed 
axial  disc  length  prescribed  by  Myrick  (1.423  in.). 

Disc  radius  was  also  incremented.  Seal  diameter  was  assumed  to 
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increase  accordingly  such  that  h,  a,  and  b  remained  constant.  Since 
Shapiro  and  Colsher's  analysis  of  Myrick's  system  also  assumed  point 
masses  at  bearing  and  disc  stations,  enlarging  the  disc  did  not  affect 
comparison  with  this  thesis'  results  provided  that  mass  values  remained 
unchanged. 

Nominal  disc  to  seal  base  clearance  (h  )  was  next  to  be 

n 

individually  varied.  The  upper  limit  for  this  parameter  was  expected  to 
be  very  restrictive  because,  beyond  a  certain  h  value,  complete  kinetic 
energy  destruction  within  the  seal  could  not  be  assured. 

Upstream  pressure  (P^)  variation  was  next.  Both  the  upper  and 
lower  limits  for  this  parameter  were  expected  to  be  sensitive.  The  upper 
limit  was  prescribed  by  choked  flow  at  seal  strip  B.  The  lower  limit  was 
restricted  by  assumed  constant  seal  exit  pressure  (Pfi) .  i.e.,  Should  PA 
become  too  small,  unrealistic  flow  reversal  would  occur.  As  with  J?  and 
h^  variations,  these  limits  were  determined  through  analysis  of  a 
complete  range  of  upstream  pressures. 

Finally,  a  and  b  seal  strip  heights  were  varied  together  to 
achieve  variable  seal  convergence  (a>  b)  and  divergence  (a<b).  For 
comparability,  average  axial  flow  area  was  maintained,  i.e..  An  increase 
in  the  "a"  dimension  was  offset  by  a  decrease  in  "b"  (Figure  3.7).  To 
establish  the  stabilizing  nature  of  a  divergent  seal,  another  speed 
variation  exercise  was  accomplished  for  an  extreme  divergent  seal 
geometry.  To  maximize  stability  threshold  increase,  PA  was  set  to 
provide  choked  flow  at  seal  strip  B .  Straight,  divergent,  and  convergent 
seal  geometry  simulations  were  also  conducted  with  exit  pressure  set  at 
one  atmosphere  (14.7  psi) .  These  simulations  indicated  whether  or  not 
this  model  is  suited  for  testing  at  higher  ambient  pressures  —  a  useful 
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feature  when  developing  high  pressure  turbine  seal  design  criteria. 


Chapter  A 


RESULTS 


A. A  Overview 

This  chapter  will  present  and  discuss  data  produced  from 
Chapter  3  simulations.  For  every  simulation,  parameter  values  will  be 
listed  followed  by  a  graphical  generalization  of  the  data  and 
observations  of  system  behavior  trends. 

Nonlinear  rotor  bearing  system  (RW5)  results  discussion  will  be 
first  and  briefest  because  of  its  limited  value  to  the  purpose  of  this 
thesis.  Linear  rotor  bearing  system  (RWEA)  results  discussions  will 
follow  and  are  limited  to  comparison  with  Shapiro  and  Colsher's  speed 
variation  results.  Isolated  seal  and  disc  linear  modeling  (RWEAP) 
results  will  be  next.  Here,  system  behavior  observations  will  take  the 
form  of  discussions  of  trends  exhibited  by  this  system  while  several 
parameters  are  individually  varied.  The  coupled  rotor  bearing  and 
labyrinth  seal  model  (RWEAS)  results  will  be  similarly  discussed.  Some 
attention  will  be  focused  upon  seal  and  rotor  bearing  systems 
interaction.  Such  discussion  provides  qualitative  understanding  only. 
Ultimately,  coupled  system  model  speed  will  be  increased  with  the  seal  at 
extreme  divergence  in  an  effort  to  prove  stability  threshold  increase. 
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Table  4.1 

Nonlinear  Rotor  Bearing  System  Parameter  Values 


Constant  Parameters 

k  =  k  =  5348.6  lbf/in. 
x  y 

o 

m  =  m  =  .0096  lbf  sec  /in 
J  L  JR 

2 

m,  =  .0875  lbf  sec  /in. 
d 

W  =  20.625  lbf 
X  =  1.603  in. 
r  =  1.081  in. 
c  =  .003  in. 

At  =  10  ^  sec 
Variable  Parameters 
Simulation  1 :  oJ  =  86  rev/sec 
Simulation  2:  u>  *=  98.3  rev/sec 
Simulation  3:  «  *  92.87  rev/sec 


=  Shaft  Stiffness 
=  Journal  Mass 
=  Disc  Mass 

=  h  Rotor  Bearing  System  Weight 
=  Journal  Length 
=  Journal  Radius 
=  Bearing  to  Journal  Radial 
Clearance 
=  Time  Step 


£  = 
i 

o 

00 

*# 

er= 

.0797 

t  = 

k 

.078, 

.07 

£  .* 

4 

.09 

H 

.074 

4.B  Rotor  Bearing  System 


4.B.I.  Nonlinear  Modeling  (RW5)  Results 


As  mentioned  in  Chapter  3,  nonlinear  rotor  bearing  system 
simulations  were  intended  only  as  a  modeling  accuracy  check.  Such 
verification  was  best  accomplished  by  operating  this  nonlinear  model  at 
or  near  Myrick's  stability  threshold. 

Simulation  1  was  conducted  below  Myrick's  stability  threshold 
of  93  rev/sec.  Damped  half  frequency  whirl  was  observed  (see  Figure 
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4.1).  Angular  velocities  for  both  journals  were  observed  to  be  identical 
for  this  and  all  other  simulations. 

Simulation  2  was  conducted  above  Myrick's  stability  threshold. 
Explosive  half  frequency  whirl  was  observed.  (see  Figure  4.2). 

Simulation  3  was  conducted  very  close  to  Myrick's  stability 
threshold.  Half  frequency  whirl  was  observed.  However,  neither  damped 
nor  explosive  whirl  could  be  confirmed  (see  Figure  4.3).  Therefore,  mass 
precession  approaching  a  limit  cycle  at  steady  state  was  assumed.  This 
condition  establishes  92.87  rev/sec  as  at,  or  very  close  to,  the  system 
stability  threshold.  Because  of  inaccuracies  inherent  in  visual  graph 
sighting,  a  limited  number  of  operating  cycles  and  computer  calculations 
this  result  was  not  expected  to  compare  exactly  with  Myrick's. 

Additionally,  lengthy  simulation  time  makes  more  accurate 
definition  of  this  stability  threshold,  using  iterative  simulations, 
unwarranted.  Therefore,  simulation  3  provided  practical  confirmation  of 
Shapiro  and  Colsher's  stability  threshold  calculation. 

4.B.2.  Linear  Modeling  (RWE4)  Results 

Linear  rotor  bearing  system  simulation  results  are  a  bit  more 
extensive  with  respect  to  nonlinear  results.  In  addition  to  a  stability 
threshold  comparison,  model  frequency  versus  rotor  speed  and  some  mode 
shapes  must  also  be  compared  with  the  results  of  Shapiro  and  Colsher. 

These  results  identify  4  primary  modes  for  Myrick's  system, 
(see  Figure  4.5a)  Mode  1  is  most  susceptible  to  self  excited  whirl.  It 
is  characterized  by  mass  displacement  due  to  both  rotor  bending  and  rigid 
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body  translation.  Mode  2  is  least  susceptible  to  any  excitation  source. 
It  is  characterized  by  rigid  body  motion  because  of  relatively  large 
journal  displacements  since  it  is  affected  predominantly  by  bearing 
coefficient  changes.  Mode  3  is  most  susceptible  to  mass  unbalance,  but 
only  marginally  affected  by  Oil  Whip  as  suggested  by  a  relatively  large 
disc  displacement.  Mode  4  is  well  damped  and  is  characterized  by  a 
combination  of  rotor  bending  and  translation. 

Rotor  speed  variation  revealed  a  stability  threshold  of  92.58 
rev/sec  (see  Figure  4.4).  This  was  indicated  by  a  near  zero  eigenvalue 
real  part  for  Mode  1.  It  does  not  differ  appreciably  from  the  nonlinear 
model  result  of  92.87  rev/sec.  Both  results  are  not  completely  accurate. 
Linear  model  accuracy  required  a  Mode  1  eigenvalue  real  part  of  exactly 
zero  and  nonlinear  accuracy  required  a  large  number  of  cycles  of 
operation.  However,  since  they  both  compare  resonably  with  Shapiro  and 
Colsher,  RW5  and  RWE4  accuracies  are  acceptable. 

Comparing  modal  frequency  vs  rotor  speed,  RWE4  results  and 
those  of  Shapiro  and  Colsher  are  in  general  agreement.  Figure  4.4a 
closely  duplicates  their  plot  (Figure  4.4b).  Per  Shapiro  and  Colsher,  a 
critical  speed  for  Modes  3  and  4  occurs  at  a  rotor  speed  of  55  rev/sec. 
RWE4  results  indicate  a  similar  critical  speed  near  55.6  rev/sec.  Both 
graphs  show  a  Mode  1  and  2  intersection  at  their  respective  stability 
thresholds  (93  and  92.58  rev/sec).  Mode  1  "growth  factors"  (eigenvalue 
real  parts)  are  also  reasonably  similar. 

Mode  shape  comparisons  are  not  as  close  as  were  previous 
comparisons  (Figure  4.5).  RWE4  rotor  speeds  cannot  be  precisely  matched 
to  those  of  Shapiro  and  Colsher  because  of  the  requirement  to  increment 
rotor  speed  by  adjusting  eccentricity  ratio  (equation  3.67).  However, 
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Figure  4.4b  -  Shapiro  &  Colsher 
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general  modal  characteristics  are  similar.  Modes  1  and  3  compare  closest 
with  Shapiro  and  Colsher.  Their  shapes  and  relative  displacements  are 
similar.  Mode  2  shapes  compare  in  that  they  both  reflect  rigid  rotor 
behavior. 

Considering  the  astonishing  similarity  between  modal  frequency 
vs  rotor  speed  plots,  these  Mode  2  differences  might  be  attributed  to 
phase  changes  occuring  between  80  and  80.13  rev/sec  for  this  mode.  Mode 
4's  shape  is  most  radically  different  from  that  of  Shapiro  and  Colsher. 
But,  since  modal  frequencies  and  growth  factors  are  similar,  these  mode 
shape  differences  might  also  be  attributed  to  phase  differences  caused  by 
non-precise  speed  duplication.  Generally,  therefore,  this  rotor  bearing 
system  model  practically  resembles  Myrick's  system  (Figure  1.8). 

4.C.  Isolated  Seal  and  Disc  (RWE4P)  -  Results 

4.C.I.  Overview 

This  section  will  discuss  results  of  parameter  variations  for 
the  isolated  and  rigidly  mounted  seal  and  disc  system.  This  system  has 
only  2  modes  because  of  such  rigid  seal  and  disc  mounting.  Consequently, 
these  modes  relate  to  flow  behavior  only.  To  identify  each  mode, 
isolating  phenomena  must  be  identified. 

When  fluid  friction  ^2h’  anc^  ^2u  coe^ic*ents  equation 
3.55)  against  disc  and  seal  surfaces  was  eliminated,  the  "Axial"  flow 
eigenvalue  was  practically  unaffected,  but,  the  "Circumferential"  flow 
eigenvalue  radically  destablized  (real  part  became  positive)  with  little 
affect  on  modal  frequency.  Since  friction  within  the  seal  affects 


Table  4.2 


Isolated  Seal  and  Disc  Nominal  Parameter  Values 


Parameter 

Program  Coding 

Value 

j e 

= 

Seal 

Strip  Separation 

CLC 

.1  in. 

r 

= 

Disc 

Radius 

RD 

10.0  in. 

PA 

s 

Seal 

Inlet  Pressure 

PSI 

(HP) 

2.55  psi 
15.00  psi 

PB 

= 

Seal 

Exit  Pressure 

PSE 

(HP) 

1.00  psi 
14.70  psi 

h 

n 

= 

Nominal  Seal  Base  to 

DC 

.200  in. 

Disc 

Clearance 

oi 

= 

Disc  Angular  Velocity 

PSDR 

579.45  rad/sec 

a 

= 

Height  of  Seal  Strip  A 

A 

.195  in. 

b 

ts 

Height  of  Seal  Strip  B 

B 

.195  in. 
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Table  4.3 

Constant  Seal  Parameter  Values 

Value 

247104  in2/sec2°R 
1.4 
560°R 

2.8  x  10  ^  lbfsec/in2 


Parameter  Program  Coding 

R  =  Gas  Constant  R 

k  =  Specific  Heat  Ratio  CK 

T  =  Seal  Fluid  Temperature  TI 

>c=  Dynamic  Viscosity  DV 


AD-AU6  774  AIR  FORCE  INST  OF  TECH  WRI GHT-PATTERSON  AFB  OH  F 

LABYKINTH  SEAL  EFFECTS  ON  ROTOR  BEARING  SYSTEM  STABILITY. (U) 
MAY  62  A  J  PAVELKO 
UNCLASSIFIED  AFIT/NR/82-10T 
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circumferential  flow  predominantly  these  modes  have  been  labled  "Axial" 
and  "Circumferential"  accordingly.  The  Axial  mode  is,  generally,  the 
more  stable  of  these  modes  as  this  section's  discussion  will  illustrate. 

The  parameter  variations  examined  in  this  analysis  are:  (1) 

seal  strip  separation  in  the  axial  direction  in  Figure  3.7),  (2)  disc 

radius  (r),  (3)  upstream  pressure  (P^),  (4)  disc  to  seal  base  nominal 

clearance  (h  ),  (5)  disc  rotational  speed  (u>) ,  and  (6)  radial  projections 
n 

of  seal  strips  toward  the  disc  (a,  b) .  Each  of  these  parameters  were 
varied  individually  leaving  all  others  at  nominal  value  (Table  4.2). 
These  nominal  values  were  selected  to  produce  realistic  results.  The 
high  pressure  simulation  (HP),  by  definition,  involved  different  nominal 
pressures. 

Nominal  values  in  Table  4.2  were  also  chosen  to  resemble  the 
seal  used  by  Wright.  The  Myrick  rotor  bearing  system  disc  radius  was 
increased  to  provide  stable  operation  with  other  nominal  parameter 
values.  Since  disc  mass  and  position  were  left  unchanged  rotor  bearing 
system  performance  was  not  affected  by  this  change. 

Throughout  this  discussion,  the  terms  "stabilize"  and 
"destabilize"  will  be  used  to  indicate  the  direction  of  change  of 
eigenvalue  position  on  a  root  locus  plot  as  a  parameter  is  incremented. 
A  "stabilizing"  trend  is  directed  toward  the  left,  "destablilizing” 
toward  the  right.  "Stable"  eigenvalues  appear  to  the  left  of  the 
imaginary  axis.  "Unstable"  eigenvalues  exist  to  the  right  of  this  axis. 

4.C.2  Seal  Strip  Separation  (/)  Variation 

Figure  4.6  presents  seal  strip  separation  variation  results. 
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For  both  Axial  and  Circumferential  modes,  stability  decreases  with 
increasing  seal  strip  separation.  This  result  appears  to  contradict  the 
expectation*  that  increasing  fluid  to  interacting  surfaces  friction  with 
increasing  seal  surface  area  would  be  stabilizing.  It  should  be  noted, 
however,  that  once  coupled  to  the  rotor  bearing  system,  this  system  will 
have  6  modes.  A  proper  stability  assessment  must  consider  all  modal 
reactions  together. 

Circumferential  modal  frequency  decreases  with  increasing  . 
This  could  represent  declining  circumferential  fluid  velocity  as  chamber 
widening  makes  this  seal  behave  more  like  a  subsonic  diffuser  by  reducing 
velocity  of  subsonic  circumferential  flow.  Corresponding  Axial  modal 
frequency  increase  suggests  that  increasing  seal  width  provides 
increasing  time  and  distance  for  acceleration  (and  deceleration)  of  axial 
flow.  For  any  of  these  simulations,  both  seal  modal  frequencies  add  to 
equal  disc  rotational  frequency  (<y) . 

Such  a  consistent  event  suggests  that  disc  rotational  frequency 
is  the  "driving"  frequency  for  this  isolated  seal  and  disc  system.  Since 
both  seal  and  disc  are  rigidly  mounted  and  seal  pressure  drop  is  constant 
only  the  Axial  and  Circumferential  "degrees  of  freedom"  are  available  to 
distribute  disc  momentum  (I^u)) .  Should  either  modal  frequency  exceed 
the  frequency  of  this  sole  momentum  source,  energy  creation  must  have 
occurred.  Obviously,  such  an  event  is  unrealistic.  Its  computer 
predicted  occurence  for  A.<  .051  in.  suggests,  possibly,  that  the 
assumption  of  complete  kinetic  energy  dissipation  within  the  seal  chamber 
is  invalid  below  this  seal  strip  separation  value. 


*  See  Chapter  1. 
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4.C.3  Disc  Radius  (r)  Variation 

Figure  4.7  presents  the  results  of  disc  radius  variation.  With 
increasing  disc  radius,  the  Axial  mode  stabilizes  while  the 
Circumferential  mode  destabilizes.  Since  Circumferential  mode 
stabilization  occurs  about  the  imaginary  axis  its  affect  on  system 
stability  is  more  important,  i.e.,  increasing  disc  radius  enhances  seal 
stability. 

As  predicted  earlier  (Chapter  1),  this  stabilization  can  be 
attributed  to  increasing  friction  between  fluid  and  adjoining  surfaces. 
Because  friction  increases  with  both  fluid  velocity  and  interacting  area 
—  both  characteristics  of  increasing  disc  radius  —  fluid  flow  through 
the  seal  should  be  stabilized. 

Another  expected  characteristic  of  increasing  r  is  increasing 
Circumferential  modal  frequency.  Increasing  circumferential  fluid 
velocity  is  intuitively  associated  with  increasing  disc  surface  speed. 
Axial  modal  frequency  must  decrease  because  of  increased  axial  flow  area 
associated  with  increasing  disc  (and  seal)  radius.  It  is  also  physically 
conceivable  that  a  gain  in  Circumferential  modal  frequency  should  offset 
a  reduction  in  Axial  modal  frequency  since  transfer  of  fluid  momentum 
involves  corresponding  changes  in  fluid  velocity.* 


*  Each  fluid  particle  has  axial  and  circumferential  velocity  components. 
Therefore,  "mass  transfer"  during  such  a  "momentum  transfer"  (P  -  mv) 
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does  not  occur. 
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4.C.4  Upstream  Pressure  (P^)  Variation 

Figure  4.8  presents  upstream  pressure  variation  results.  Both 

Axial  and  Circumferential  modes  destabilize  when  P.  is  increased. 

A 

Increasing  upstream  pressure  implies  increasing  force  applied  against 
seal  and  disc.  Since  both  disc  and  seal  are  physically  constrained  only 
flow  modes  can  react  to  this  increasing  force.  Such  a  reaction  should 
reflect  reduced  stability. 

Unfortunately,  actual  destabilization  (left  to  right  crossing 
of  imaginary  axis)  is  not  illustrated  here.  At  P^  =  2.55  psi,  flow 
through  seal  strip  B  (Figure  3.7)  chokes.  Further  increase  in  PA>  until 
flow  through  seal  strip  A  chokes,  should  destabilize  the  Circumferential 
Mode.  In  addition  to  this  upper  realistic  limit,  a  lower  realistic  limit 
was  also  encountered. 

Circumferential  modal  frequency  exceeds  disc  rotational 
frequency  when  PA  <  1.3108  psi.  As  explained  in  4.C.2,  such  a  condition 
is  unrealistic.  This  condition,  possibly,  corresponds  to  a  reversal  of 
flow  direction  thereby  violating  a  system  constraint  —  PA  >  Pp. 

4.C.5.  Disc  to  Seal  Base  Nominal  Clearance  (h  )  Variation 

n 

Figure  4.9  presents  results  of  h^  variation  simulations. 
Within  realistic  limits,  both  Axial  and  Circumferential  modes  are 
stabilized  by  increasing  h^.  This  result  might  be  attributed  to  a 
reduced  rate  of  change  of  chamber  pressure.  As  hR  increases,  axial  fluid 
velocity  and  modal  frequency  decrease  th  ereby  reducing  rate  of  change  of 


chamber  pressure  (equation  3.51  -  £P  coeff icient) . 

This  effect  is  valid  between  narrow  limits  only.  The  obvious 
lower  limi'  is  h^  =  a  and/or  h^  =  b  resulting  in  complete  cessation  of 
flow.  The  upper  limit  is  h^ St  .27  in.  (interpolated  result).  This  limit 
corresponds  to  an  unrealistic  Circumferential  modal  frequency  greater 
than  u)  .  Consequently,  this  limit  corresponds  to  a  clearance  so  great 
chat  the  assumption  of  complete  kinetic  energy  destruction  within  the 
seal  chamber  is  violated. 

4.C.6.  Disc  Angular  Velocity  (ia)  Variation 

Figure  4.10  presents  disc  angular  velocity  variation  results. 
Both  Axial  and  Circumferential  modes  stabilize  with  increasing^.  This 
stabilization  can,  again,  be  attributed  to  increasing  friction  between 
seal  fluid  and  interacting  surfaces  with  increasing  disc  surface  speed. 

4.C.7.  Seal  Divergence  (a >  b)  Variation 

Figure  4.11  presents  seal  divergence  variation  results.  Both 
seal  modes  are  stabilized  by  increasing  seal  divergence.  Decreasing 
nominal  chamber  pressure  with  increasing  seal  divergence  suggests  that 
greater  exit  with  respect  to  inlet  mass  flow  area,  reducing  chamber 
pressure,  is  responsible  for  this  stabilizing  trend. 

4.C.8  Seal  Convergence  (a  <  b)  Variation 

Figure  4.12  presents  seal  convergence  variation  results.  Seal 
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convergence  destabilizes  both  seal  modes.  Convergence  results  are 
expected  to  be  equal  and  opposite  to  seal  divergence  results. 
Unfortunately,  identical  support  pressures  (P  and  P  )  for  both 

A  Sj 

simulations  causes  choked  flow  at  seal  strip  B  and,  therefore,  constant 
nominal  chamber  pressure  with  increasing  seal  convergence.  However, 
evident  reduced  stabilization,  with  respect  to  seal  divergence,  can  be 
attributed  to  rate  of  change  of  chamber  pressure  per  equations  3.48. 
These  mass  flow  rate  variation  relations  are  affected  by  axial  flow  area. 
They  are  also  strong  coefficients  of  S P  in  continuity  equation  3.51. 

4.C.9  Variable  Seal  Geometry  at  High  Chamber  Pressure 

The  results  in  Figure  4.13  indicate  that  straight  seal  geometry 
is  not  as  stable  as  convergent  or  divergent  geometries  for  seal  modes  at 
high  chamber  pressure.  This  result  might  indicate  that  any  unequal  seal 
geometry  is  stabilizing  at  high  pressure  and/or  low  axial  fluid  velocity. 
However,  it  is  believed  that,  for  flow  Mach  Numbers  closer  to  1.0,  these 
results  would  be  similar  to  those  in  Figures  4.11  and  4.12. 

4.C.10  General  Trends 

The  most  obvious  trend  among  these  RWE4P  simulations  is  the 
Circumferential  mode's  role  as  the  least  stable  seal  mode.  In  every 
simulation,  it  is  either  the  least  stable  or  the  only  unstable  seal  mode. 
Eigenvector  analysis  further  supports  this  observation.  For  all 
simulations,  the  largest  eigenvector  component  is  the  "A^"  component 
(equation  3.56b).  This  cosine  related  circumferential  velocity 
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coefficient  exceeds  axial  pressure  components  (A^  and  in  equation 

4 

3.56a)  by,  at  least,  a  factor  of  10  .  Consequently,  the  Circumferential 
mode,  as  the  most  volatile  seal  mode,  deserves  maximum  attention  when 
stability  is  discussed.  Additionally,  its  modal  frequency  always  exceeds 
that  of  the  Axial  mode. 

One  more  pattern  common  to  all  RWE4P  simulations  is  that  of 
Axial  and  Circumferential  modal  frequencies  adding  to  equal  disc 
rotational  frequency  (u)  ) .  Explained  earlier,  this  criterion  is 
responsible  for  locating  many  realistic  limits  for  seal  parameters. 

Of  all  parameter  variations  performed  with  this  isolated  disc 
and  seal  system,  increasing  seal  divergence  and  disc  speed  are  the  most 
stabilizing.  Increasing  disc  radius  also  has  an  overall  stabilizing 
affect  on  the  seal.  The  stabilizing  effects  of  increasing  both  to  and  r 
are  due  to  consequent  increasing  friction  between  seal  fluid  and 
interacting  surfaces.  Although  seal  base  to  disc  nominal  clearance 
increase  also  stabilizes  both  seal  modes,  its  limited  realistic  range 
made  its  usefullness  as  a  stabilizing  design  parameter  unpromising. 

Because  most  RWE4P  results  have  satisfactory  physical 
explanations  they  can  serve  somewhat  as  a  modeling  accuracy  verification. 
Their  primary  purpose,  however,  is,  still,  to  provide  information 
describing  affects  that  individual  seal  parameters  have  on  seal 
stability.  In  the  next  section,  these  parameter  variations  will  be 
evaluated  in  the  context  of  a  coupled  seal  and  rotor  bearing  system. 
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4.D  Coupled  Rotor  Bearing  and  Seal  System  (RWE4S)  -  Results 

4.D.1  Overview 

Now  that  both  rotor  bearing  and  labyrinth  seal  models  have  been 
operationally  proven,  coupled  system  reactions  to  seal  parameter 
variations  can  be  discussed.  First,  the  behavior  of  Myrick's  4  modes 
will  be  presented  graphically.  All  modes  are  identified  by  comparing 
RWE4S  eigenvalues  with  RWE4  and  RWE4P  eigenvalues  obtained  under  the  same 
operating  conditions.*  When  an  RWE4S  eigenvalue  pair  more  closely 
equals  one  RWE4  or  RWE4P  eigenvalue  pair  than  any  other,  it  is  assigned 
that  RWE4  or  RWE4P  mode's  designator.  Secondly,  differences  in  behavior 
for  these  modes  between  coupled  and  rotor  bearing  (RWE4)  systems  will  be 
discussed.  Such  a  discussion  will  explain  seal  influence  on  the  rotor 
bearing  system.  Finally,  rotor  bearing  system  influence  on  seal  modes 
will  be  discussed  by  relating  coupled  system  results  with  those  for  the 
isolated  seal  and  disc. 

The  format  of  this  discussion  will  parallel  that  used  for 
discussing  RWE4P  results.  The  same  seal  and  disc  parameters  will  be 
varied  individually.  Corresponding  plots  for  Modes  1  through  4  will  be 
presented  and  explained. 

4.D.2  Seal  Strip  Separation  ( k )  Variation 

Figure  4.14  presents  seal  strip  separation  variation  results 

*  Parameter  values  in  Tables  4 . 1  and  4 . 2 


for  Myrick's  first  4  modes. 

Mode  1  behaves  as  expected  for  &  values  that  are  near  nominal 
(Table  4.2).  Increasing  H  beyond  .1  in.  (nominal)  stabilizes  this 
fundamental  Myrick  mode.  This  event  might  be  interpreted  as  a  "stability 
transfer"  between  Mode  1  and  destabilizing  seal  modes  (Figure  4.6). 
However,  Mode  1  promptly  begins  to  destabilize  for  Jl  >  .1255  in., 
indicating  an  upper  limit  for  such  a  stability  transfer.  This  stability 
turning  point  also  corresponds  closely  to  actual  destabilization  of  the 
Circumferential  seal  mode.  Another  eigenvalue  real  part  direction  change 
toward  increasing  stability  at  Jl  =  .2002  in.  could  be  either  an  aberation 
or  an  indication  that  seal  strip  separation  growth  cannot,  by  itself, 
destabilize  Mode  1.  However,  an  absence  of  eigenvalue  real  part 
direction  change  for  seal  modes  at  *  .2002  in.  suggests  that  Mode  1 
destabilization  will  continue  with  continued  increase  in  Jl  .  Further 
speculation  of  such  a  stability  transfer  must  await  presentation  of 
remaining  J?  variation  results. 

The  remarkable  stability  of  Mode  2  justifies  Shapiro  and 
Colsher's  high  appraisal  of  this  mode  (Figure  4.5).  Despite 
computational  accuracy  to  7  digits  (single  precision),  no  change  in 
stability  or  modal  frequency  is  noticeable. 

Mode  3's  behavior  is  similar  to  that  of  Mode  1,  but,  with 
larger  real  part  variation  and  reversal  of  modal  frequency  change 
direction.  It,  too,  stabilizes  with/?  increase  up  to  .1255  in.  However, 
its  modal  frequency  decreases  for  X  >  .1006  in.,  after  a  gradual  increase 
for  X  <  .1006  in.,  and  then  begins  another  gradual  increase  for  X  >  .1504 
in. 

Mode  4  destabilizes  with  increasing  Jl .  Its  range  of  eigenvalue 
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real  part  variation  is  between  that  for  Modes  1  and  3.  At  J?  =  .1255  in., 
its  rate  of  destabilization  decreases.  Its  modal  frequency  also  begins 
to  decrease  at  this  point  after  a  previous  gradual  increase. 

For  rotor  bearing  system  Modes  1,  3,  and  4,  .1255  in. 
represents  a  transition  point.  Modes  1  and  3  begin  destabilization  while 
Mode  4  continues  its  destabilization.  This  point  also  roughly  corresponds 
to  actual  Circumferential  seal  mode  destabilization.  It  is  plausible, 
therefore,  that  the  straight  geometry  seal  has  a  stabilizing  effect  on 
the  rotor  bearing  system's  most  delicate  modes  (1  and  3)  until  the  seal 
itself  becomes  unstable.  Such  an  effect  might  also  be  described  as  a 
"stability  transfer"  from  seal  modes  to  Modes  1  and  3.  This  suggests 
that  unstable  circumferential  flow  within  the  seal  acts,  through 
friction,  to  destabilize  all  rotor  bearing  system  modes.  Therefore, 

1255  in.  represents  an  upper  design  limit  for  RWE4S  straight  seal 
stability. 


4.D.3  Disc  Radius  (r)  Variation 

Figure  4.15  presents  r  variation  results  for  Modes  1  through  4. 

Mode  1  stabilizes  normally  with  increasing  r,  i.e.,  it  follows 
a  smooth  trend  without  observed  aberrations.  Its  modal  frequency  also 
increases  with  increasing  r. 

Mode  2  is,  again,  unaffected.  However,  Mode  3  behaves  almost 
in  an  equal  and  opposite  manner  with  respect  to  JL  variation  results. 
Here,  r  ■  10  in.  (nominal  value)  is  close  to  a  transition  toward 
destabilization.  Since  this  value  also  represents  near  instability  for 
the  seal's  Circumferential  mode  (Figure  4.7)  a  "stability  transfer" 
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theory  similar  to  that  proposed  for  X  variation,  again,  appears 
plausible . 

Mode  4  destabilizes  with  increasing  r.  Its  modal  frequency 
simultaneously  increase  ?  The  value  r  =  10  in.  appears  to  have  no 
special  significance  for  this  mode.  However,  this  root  locus  plot  is 
more  erratic  than  that  for  Mode  1,  for  either  seal  mode,  or  for  Mode  4 
during  A.  variation.  Considering  also  Mode  1's  relatively  smooth 
stabilizing  behavior  and  the  similarity  in  form  between  Mode  3  in  both  r 
and  X  variations,  increasing  disc  radius  seems  to  have  an  approximately 
opposite  effect  on  rotor  bearing  system  behavior  than  does  increasing 
seal  strip  separation. 

Since  increasing  both  A  and  r  should  increase  disc  to  seal 
fluid  friction  and  since  the  rotor  bearing  system  disc  is  modeled  as  a 
point  mass,  such  friction  is  not  responsible  for  this  difference. 
Instead,  some  other  feature  of  the  Circumferential  seal  mode  should  be 
considered  responsible.  This  is  the  only  seal  mode  that  destabilizes 
with  increasing  A  while  stabilizing  with  increasing  r. 

4.D.4  Upstream  Pressure  (PA>  Variation 

Figure  4.16  presents  PA  variation  results  for  Myrick's  modes. 
These  results  are  much  less  erratic  than  those  for  A  anc*  r  variations. 

Mode  1  stabilizes  with  increasing  PA-  Its  modal  frequency 
simultaneously  decreases  slightly.  Mode  2  is  unaffected.  Mode  3 
destabilizes  with  increasing  P^.  Mode  4  stabilizes  with  increasing  PA 
while  its  modal  frequency  decreases. 

Increasing  pressure  difference  across  the  seal  is  intuitively 
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believed  to  be  destabilizing.  Accordingly,  Mode  1  should  be  destabilized 
by  such  increasing  pressure  since  it  is  the  most  susceptible  mode  to  self 
excited  instability.  However,  these  results  challenge  such  a  belief. 
This  surprising  result  is  significant  in  that  increasing  upstream 
pressure  might  help  reduce  the  nemesis  of  Mode  1  —  speed  induced  or 
self  excited  instability.  By  increasing  upstream  pressure  when  rotor 
speed  approaches  the  stability  threshold.  Mode  1  instability  might  be 
postponed. 


4.D.5  Seal  Base  to  Disc  Nominal  Clearance  (h  )  Variation 

n 

Results . 

Figure  4.17  presents  hR  variation  results.  As  in  isolated  disc 

and  seal  simulations,  an  upper  realistic  limit  is  encountered.  Root 

locus  plots  reflect  unrealistic  eigenvalues  with  dotted  lines. 

Mode  1  is  stabilized  by  increasing  h^  within  realistic  limits 

while  its  modal  frequency  is  decreased.  Mode  2  is  finally  affected  by  a 

seal  parameter  variation.  It  is  slightly  stabilized.  Its  modal 

frequency  is  increased  with  increasing  h^.  Mode  3  is  destabilized  with 

increasing  h  while  its  modal  frequency  is  increased.  Mode  4  is  also 
n 

destabilized.  Its  modal  frequency  is  decreased. 

4.D.6  Disc  Angular  Velocity  («))  Variation 

Figure  4.18  presents  disc  speed  variation  results.  Mode  1 
demonstrates  a  stability  threshold  (92.69  rev/sec)  slightly  higher  than 
the  RWE4  result  (92.56  rev/sec).  It  also  exhibits  a  stabilizing  trend 
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for co  values  less  than  38.22  rev/sec  (240  rad/sec).  Its  modal  frequency 
increases  smoothly  with  increasing^. 

Mode  2  behaves  similarly.  Its  transition  toward 
destabilization  with  increasing  begins  at  about  69  rev/sec  (435 
rad/sec) .  Although  beyond  range  of  physical  operation,  its  stability 
threshold  is  about  464  rev/sec  (2915  rad/sec).* 

Mode  3  behaves  erratically  during  speed  variation.  Initially, 
it  stabilizes  with  u)  increases  up  to  76  rev/sec  (477  rad/sec).  This 
eigenvalue  plot  changes  direction  3  more  times  before  starting  a 
destabilizing  journey  that,  apparently,  never  quite  reaches  the  imaginary 
axis  (eigenvalue  real  part  at  1392  rev/sec  =  -.63,  corresponding  to  E.  = 
.005).  Such  erratic  behavior  supports  Mode  3's  reputation  as  a  lightly 
damped  mode  readily  destabilized  by  mass  unbalance  but  not  by  increasing 
speed , 

Mode  4  is  stabilized  by  increasing^  during  the  useful  range  of 
this  simulation.  However,  like  Modes  1  and  2,  it  also  turns  toward  the 
imaginary  axis.  This  turning  point  occurs  at  about  126  rev/sec  (790 
rad/sec) .  Its  fictitious  stability  threshold  is  just  below  that  of  Mode 
2  thereby  supporting  eventual  convergence  of  Modes  2  and  4  with 
increasing u)  (Figure  4.4  —  RWE4  Results).  Similarly,  at  another  large  w 
value,  Modes  1  and  3  converge  —  Mode  1,  presumably,  transverses  the 
imaginary  axis  several  times  while  Mode  3  approaches  it  asymptotically. 


*  Mode  1  periodically  becomes  unstable  and  then  stable  as  speed  increases 
toward  this  "threshold".  Therefore,  physical  failure  is  expected  well 
before  this  speed  is  achieved. 


4.D.7  Seal  Divergence  Variation 


Figure  4.19  presents  seal  divergence  results.  As  noted, 
simulation  1  represents  a  straight  geometry  seal.  Simulations  2  and  3 
represent  divergent  seal  geometries.  Total  flow  area  through  both  seal 
strips  is  maintained. 

Mode  1  is  stabilized  by  increasing  seal  divergence.  Its  modal 

frequency  is  also  increased.  Mode  2,  again,  proves  itself  resistant  to 

change.  Mode  3  reexhibits  its  penchant  for  erratic  behavior.  It  is 

stabilized  by  moderate  seal  divergence,  but,  destabilized  by  continued 

divergence.  Mode  4  is  destabilized  by  increasing  seal  divergence. 

Since  stabilization  of  Mode  1  is  required  to  increase  stability 

threshold  seal  divergence  is,  apparently,  beneficial.  Since  degree  of 

* 

seal  divergence  is  physically  limited  ,  actual  destabilization  of  Mode  3 
or  4  is  unlikely  with  increasing  seal  divergence. 

4.D.8  Seal  Convergence  Variation 

Figure  4.20  presents  seal  convergence  variation  results. 
Conditions  stated  for  seal  divergence  apply  here  also. 

As  expected,  Mode  1  is  destabilized  by  increasing  seal 
convergence.  Its  small  degree  of  dtstabilization  relative  to  variable 
divergence  induced  stabilization  (Figure  4.19)  can  be  attributed  to 
constant  nominal  chamber  pressure  caused  by  choked  flow  at  seal  strip  B. 


*  .  19  in.  <.  a,b  <  h 
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Mode  3  is  stabilized  by  increasing  seal  convergence  while  mode 
4  is  destabilized.  Again,  rates  of  eigenvalue  real  part  change  relative 
to  variable  seal  divergence  simulations  are  marginal  due  to  constant 
nominal  chamber  pressure. 

4.D.9  Variable  Seal  Geometry  at  High  Chamber  Pressure 


Figure  4.21  presents  high  chamber  pressure  results.  Here,  seal 
geometry  varies  from  convergent  to  divergent  while  relative  pressure  drop 
across  the  seal*  is  maintained  at  a  relatively  low  level.  Nominal 
chamber  pressure  is  increased  approximately  by  a  factor  of  10.  The 
stability  reversals  evident  for  seal  modes  (Figure  4.13)  are  not 
exhibited  by  rotor  bearing  system  modes. 

Mode  1  is  stabilized  by  decreasing  seal  convergence  (or 
increasing  divergence)  thereby  maintaining  predicted  behavior  (Figure 
4.19).  Mode  3  also  stabilizes  as  expected.  Mode  4  satisfies  the 
prediction  of  Figure  4.19  by  destabilizing  with  increasing  divergence. 

It  should  be  remembered  that  these  results,  particularly  Mode 
l's,  are  of  primary  interest  when  considering  stability  threshold 
enhancement.  Here,  predicted  Mode  1  stabilization  with  diverging  seal 
geometry  is  encouraging.  Unfortunately,  the  seal  Circumferential  mode  is 
unstable  during  all  of  these  high  pressure  simulations  making 
experimental  verification  impossible.  Additionally,  more  simulations  at 
other  pressure  (PA»  Pg)  and  pressure  drop  (P  -Pg)  values  are  needed  to 
provide  confidence  in  maintaining  this  assertion  of  improved  stability  at 
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high  chamber  pressure  using  divergent  seal  geometry. 

4.D.10  Speed  (ti>)  Variation  at  Extreme  Seal  Divergence 

Figure  4.22  presents  root  locus  results  for  Myrick  modes  at 
extreme  seal  divergence  and  for  variable  rotor  speeds  starting  near  the 
Mode  1  stability  threshold  of  Figure  4.18.  This  simulation  is  intended 
to  determine  a  revised  stability  threshold  due  to  seal  divergence  and  to 
verify  improved,  divergence  induced,  coupled  system  stability  per  Figure 
4.19. 

Mode  l's  stability  threshold  is  increased  to  93.13  rev/ sec 
(585.16  rad/sec)  because  of  seal  divergence.  This  represents  a  .47% 
increase  from  the  revised  stability  threshold  achieved  with  a  straight 
seal  (Figure  4.18)  and  a  .6%  increase  relative  to  the  stability  threshold 
for  Myrick's  rotor  bearing  system  without  seal  (Figure  4.4). 

Stability  of  Modes  2  and  3  is  marginally  degraded,  but  not 
enough  to  discourage  further  stabilization  of  Mode  1.  Mode  4  continues 
its  stabilizing  behavior  with  increasing (Figure  4.18). 

4.D.11  Seal  and  Rotor  Bearing  System  Interactions 

4.D.11.A  Seal  Influence  on  Rotor  Bearing  System  Modes 

Myrick  modes  are  not  radically  affected  by  coupling  with  a 
labyrinth  seal.  Generally,  points  on  root  locus  plots  in  Figures  4.14 
through  4.22  are  affected  equally  by  seal  modes.  Although  such  affects 
cannot  be  positively  attributed  to  any  particular  seal  mode,  the  relative 
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volatility  of  the  Circumferential  mode  suggests  its  predominance. 
Therefore,  curve  shapes  and  relative  distances  among  points  remain 
practically  unchanged. 

Mode  1  is  stabilized  by  the  seal.  The  seal  also  acts  to 

increases  its  modal  frequency. 

As  most  RWE4S  simulations  reveal,  the  seal  has  little 
perceptible  effect  on  Mode  2.  Only  speed  variations  (not  a  seal 
parameter)  significantly  affects  Mode  2  (Figures  4.18  and  4.22). 

Mode  3  is  destabilized  and  has  decreasing  modal  frequency  for 
W^  477  rad/sec  (76  rev/ sec)  due  to  seal  influence.  For  w>  477  rad/sec, 
Mode  3  is,  conversely,  stabilized  and  has  increasing  modal  frequency  as  a 
result  of  seal  influence. 

Mode  4  is  destabilized  by  the  seal.  Its  modal  frequency  is 

increased. 


4.D.11.B  Rotor  Bearing  System  Influence  on  Seal  Modes 

Rotor  bearing  system  influence  on  seal  modes  is  similar  to  seal 
influence  on  rotor  bearing  system  modes  in  that  RWE4P  seal  mode  curve 
shapes  and  relative  distances  among  points  are  unaffected.  The  Axial 
mode  is  stabilized  while  the  Circumferential  mode  is  destabilized  when 
seal  and  rotor  bearing  systems  are  coupled.  Again,  because  rotor  bearing 
system  modes  were  not  isolated,  no  particular  mode  can  be  identified  as 
having  a  predominant  affect  on  these  seal  behavior  changes. 

Another  such  seal  behavior  change  deals  with  affects  on  seal 
modal  frequencies.  When  seal  modes  are  stable,  the  rotor  bearing  system 
acts  to  reduce  modal  frequencies  for  both  seal  modes.  Such  an  event 
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suggests  "frequency  transfer"  from  the  seal  to  the  rotor  bearing  system 
since  seal  modal  frequencies  now  add  to  be  less  thanu)  .  However,  when 
the  seal  (Circumferential  mode)  is  unstable,  seal  modal  frequencies  add 
to  exceed  u> .  In  this  case,  Axial  (still  stable)  modal  frequency  is  still 
reduced,  but  Circumferential  modal  frequency  increases.  Since  this 
result  was  not  encountered  during  isolated  seal  and  disc  simulations  this 
excess  frequency  must  be  assumed  to  have  the  form  of  disc  whirl.  Because 
both  systems  are  coupled  with  an  unstable  seal  mode  this  whirl  must  be 
explosive — a  dangerous  but  still  realistic  condition! 

4 . E  Summary 

RW5  and  RWE4  rotor  bearing  system  results  compare  favorably 
with  those  of  Myrick  and  Shapiro  and  Colsher.  Such  success  establishes 
Chapter  3's  analytical  rotor  bearing  system  model  as  both  accurate  and 
similar  to  Myrick' s  experimental  model. 

RWE4P  results  reveal  increasing  disc  radius,  disc  rotational 
speed,  seal  base  to  disc  nominal  clearance  and  seal  divergence  as  having 
stabilizing  influences  on  seal  fluid  flow.  Conversely,  increasing  seal 
strip  separation,  chamber  pressure,  and  seal  convergence  destablize  fluid 
flow  through  the  seal.  At  high  pressure  and  low  axial  flow  rate,  both 
diverging  and  converging  seal  geometries  stabilize  seal  modes  relative  to 
straight  seal  geometry.  These  results,  generally,  have  plausible 
physical  explanations  thereby  verifying  seal  modeling  accuracy. 

Coupled  system  results  also  compare  favorably  with  earlier 
predictions.  Because  Mode  1  is  most  susceptible  to  self  excited 
instability,  parameter  variation  affects  on  it  determine  coupled  system 
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stability. 

Mode  1  is  stabilized  by  increasing  disc  radius,  seal  pressure 
drop,  nominal  seal  base  to  disc  clearance,  and  seal  divergence.  It  is 
destabilized  by  increasing  seal  convergence.  In  the  vicinity  of  nominal 
parameter  values  (Table  4.2),  increasing  seal  strip  separation  stabilizes 
Mode  1.  However,  an  unstable  seal,  with  further  increase  in  ^  ,  soon 
causes  Mode  1  to  destabilize.  The  destabilizing  effect  of  increasing 
disc  speed  on  Mode  1  is  caused  primarily  by  Oil  Whip.  Straight  seal 
geometry  postpones  the  onset  of  Oil  Whip  marginally. 

Similarly,  diverging  seal  geometry  increases  system  stability 
threshold.  Although  still  marginal  in  degree,  this  increase  satisfies 
expectations  in  that,  amount  of  stability  threshold  increase  exceeds  that 
of  both  straight  and  convergent  seals.  This  marginal  increase  cannot  be 
considered  within  range  of  modeling  error,  indicated  in  ^.B.2,  because 
RWE4  and  RWE4S  both  contain  the  same  rotor  bearing  system  programming. 
Additionally,  by  increasing  stability  threshold  the  least  relative  to 
seal  absence,  convergent  seal  geometry  results  futher  reinforce  this 
prediction  of  increasing  stability  threshold  with  increasing  seal 
divergence.  This  trend  appears  to  apply  also  at  increased  chamber 
pressure  and  reduced  axial  flow  rate.  Unfortunately,  additional  data  at 
different  combinations  of  such  higher  chamber  pressure  and  lower  flow 
rates  are  needed  for  verification.  Such  verification  was  not 
accomplished  using  this  model  because  of  an  unstable  Circumferential  seal 


mode  at  high  chamber  pressure. 


Chapter  5 


CONCLUSIONS 


The  principle  aim  of  this  thesis  has  been  achieved.  Rotor 
bearing  system  stability  has  been  improved  by  divergent  seal  geometry. 
Although  such  improvement  was  marginal  (.47%),  it  should  be  noted  that  a 
specific  maximized  stability  threshold  increase  was  not  sought.  System 
adjustment  to  achieve  such  a  threshold  is  warranted  only  when  dealing  with 
models  directly  proportionate  with  operational  steam  turbines. 

Seal  parameter  affects  on  system  stability  were  discovered  and 
presented  for  possible  later  use  in  developing  seal  design  criteria. 
Such  information  might  be  used  to  optimize  stability  of  both  seal  and 
rotor  bearing  system  modes  toward  the  goal  of  turbine  efficiency 
enhancement. 

These  specific  and  notable  conclusions  were  developed: 

1)  Performances  of  nonlinear  (RW5)  and  linear  (RWE4)  analytical  rotor 
bearing  system  models  are  similar  to  that  of  Myrick's  experimental  model. 

2)  Labyrinth  seal  stability  is  governed  primarily  by  its  Circumferential 
mode. 

3)  Straight  seal  geometry  increases  rotor  bearing  system  stability. 

*4)  Divergent  seal  geometry  increases  rotor  bearing  system  stability 
threshold  beyond  that  of  straight  seal  geometry. 

5)  Convergent  seal  geometry  destabilizes  the  Myrick  rotor  bearing  system 
relative  to  straight  and  divergent  seal  geometries. 

6)  This  coupled  rotor  bearing  and  labyrinth  seal  system  model  is  not 
suited  for  experimental  testing  at  atmospheric  pressure. 


APPENDIX  A 


A. 1  Nonlinear  Rotor  Bearing  System  Programming 


PROGRAM  RW5 (RW5 .OUTPUT, TAPE5=RW5 , TAPE6=0UTPUT) 

DIMENSION  X(20) , F(20) ,TS(100) ,X3(100) ,K(7) ,X1(100) ,X2(100) 
COMMON/DR/SKX, SKY , BC , W, PI , SPI , DVR , DVL , CL , BR, PSDR, PER, PEL 
1 , ERR, ERL, XS (20) , H , DM, BM, FXLS , FYLS , FXRS , FYRS 
DATA  X,F,TS ,X5,XS/452*0. /  $DATA  K/7*0/ 

C  THIS  PROGRAM  SIMLLATES  ROTOR  BEARING  SYSTEM  OPERATION  USING  FINITE 
C  DIFFERENCE  TIME  INTEGRATION 

C  H-TIME  STEP=l.E-05  SEC 

C  DV= DYNAMIC  VISCOSITY  OF  BEARING  FLUID=5 ■ E-06LBFSEC/IN2 

C  CL3 JOURNAL  LENGTH=1 .06 3 IN 

C  W»Sj  (  JOURNAL+D  IS  C  )  WEIGHT 

C  BC= RADIAL  BEARING  CLEARANCE3 . 003IN 

C  BR= JOURNAL  RADIUS3 1 . 081 IN 

C  SKX=X  DIRECTION  ROTOR  SHAFT  STIFFNESS=5348. 6LBF/IN 

C  SKY=Y  DIRECTION  ROTOR  SHAFT  STIFFNESS=5348 . 6LBF/IN 

C  DM3DISC  MASS3 .0875LBFSEC2/IN 

C  BM3 JOURNAL  MASS=.0096LBFSEC2/IN 

C  ERI=INITIAL  JOURNAL  ECCENTRICITY  RATIO 

C  ERS=STEADY  JOURNAL  ECCENTRICITY  RATIO 

C  PEI3INITIAL  ATTITUDE  ANGLE 

C  IH=TIME  STEP  COUNTER 

C  ITP=PRINT  COUNTER 

C  MT=MAX  NUMBER  OF  TIME  STEPS 

C  NN-NUMBER  OF  STATE  VARIABLES3 12 

C  M3NUMBER  OF  PLOT  POINTS 

C  NP=PLOTTING  OPTION (MINIPLT) 

READ (5,1)H,DV,CL,W,BC,BR,SKX,SKY, DM , BM , ERI , PEI , ERS 

READ(5,2)IH,ITP,MT,NN,M,NP 

PI=3. 1415927 

SPI=9. 8696044 

C  INITIALIZE  ECCENTRICITY  RATIO 

XS(15)=ERS  $  XS(13)=ERS 
DVR=5 . DE-06  $DVL=5 . DE-D6 
Dl*(l ,-(XS(15)**2) ) 

D2=XS(15) / (Dl**2) 

D3=((16./SPI)-1.) 

D4= (D3* (XS ( 1 5 ) **2) +1 . ) 

C  CALCULATE  JOURNAL  ANGULAR  VELOCITY 

PSDR34.*W*(BC**2)*(D1**2)/(PI*DVR*CL*XS(15)*SQRT(D4)*BR) 

C  CALCULATE  ATTITUDE  ANGLE 

XS(16)=ATAN((PI/4.)*SQRT(D1)/XS(15)) 

C  LEFT  BEARING 

21  D1-(1.-(XS(13)**2)) 

D23XS(1 3) / (Dl**2) 

D3=((16./SPI)-l.) 

D4=(D3*(XS(13)**2)+1.) 

PSDL=4.*VJ*(BC**2)*(D1**2)/(PI*DVL*CL*XS(13)*SQRT(D4)*BR) 


150 


151 


XS(14)=ATAN( (PI/4 . )*SQRT(Dl) /XS (13) ) 
T=0.0 

STEADY  CONDITIONS 
XS(3)=XS(13)*BC*SIN(XS(14) ) 

XS (4) =-XS ( 1 3) *BC*COS (XS (14) ) 
XS(5)=XS(3) 

XS(6)=XS(4) 

XS(1)=XS(3) 

XS(2)=XS(4) 

STEADY  BEARING  FILM  FORCES 
FXLS=W*SIN(XS(14) ) 

FYLS=W*COS(XS(14) ) 

FXRS=W*SIN(XS(14)) 

FYRS=W*COS(XS(14)) 
WRITE(6,6)T,H,XS(15) ,XS(16) ,PSDR 
WRITE(6?8)XS(13) ,XS(14) ,PSDL 
WRITE (6 ,7)(XS(I),I=1 ,NN) 

INITIAL  CONDITIONS 

X(3)=ERI*BC*SIN(PEI) 

X(4)=-ERI*BC*COS(PEI) 

X(5)=X(3) 

X(6)=X(4) 

X(1)=X(3) 

X(2)=X(4) 

X(13)=SQRT( (X(3)**2)+(X(4)**2))/BC 
X(15)=SQRT( (X(5) **2)+(X(6) **2) ) /BC 
X(14)=PEI 
X(16)*PEI 

WRITE(6,6)T,H,X(15) ,X(16) ,PSDR 
WRITE(6,8)X(13) ,X(14) ,PSDL 
WRITE(6,7)(X(I),I=1,NN) 

KI=0 

DO  10  11*1 ,MT,ITP 
KI*KI+1 
TS(KI)=T 
X1(KI)=X(1) 

X2(KI)=X(2) 

X3(KI)=X(3) 

CALCULATE  STATE  VARIABLE  VALUES 
DO  11  JJ*1,ITP,IH 
CALL  RK4(F,X,T,NN,H) 

11  CONTINUE 

WRITE(6,6)T,H,X(15) ,X(16) ,PSDR 
WRITE (6 ,8)X(1 3) ,X(14) ,PSDL 
WRITE (6 , 7) (X(I) ,1=1 ,NN) 

10  CONTINUE 
33  K(l)=l  $  K(3)=NP 

K(6)=10HX1  VS  T 
PLOT  STATE  VARIABLE  VALUES 
CALL  MINIPLT(TS,X1 ,M,K) 

K(6)=10HX2  VS  T 
CALL  MINIPLT(TS ,X2 ,M,K, ) 

K(6)*10HX2  VS  XI 
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CALL  MINIPLT(X1 ,X2,M,K) 

K(6)=10HX3  VS  T 
CALL  MINIPLT (TS,X3,M,K) 

7  FORMAT (  E23.6) 

6  FORMAT (  / , 5X, 'T  =  ’ ,E10.3,5X, *H  =  ' ,E15.6,3X, 'ERR=' ,E12.6,3X, 
l'PER=',E12.5,3X, ' PSDR= ' ,E12 . 5) 

S  FORMAT(46X, "ERL= ' .E12.6.3X, 'PEL-' ,E12.5,3X, ' PSDL= ' ,E12. 5, /) 

1  FORMAT (E20. 7) 

2  FORMAT(IIO) 

STOP 

END 

SUBROUTINE  DERIV(F,X,T) 

DIMENSION  X(20) ,F(20) 

COMMON/DR/ SKX , DKY , BC , W , PI , SPI , DVR  ,  DVL , CL , BR , PSDR , PER, PEL 
1,ERR,ERL,XS(20),H,DM,BM,FXLS,FYLS,FXRS,FYRS 
C  ECCENTRICITY  AND  ATTITUDE  ANGLE  INCREMENTS 
E=X(I3) *BC  $  ES=E**2 

F(13)=(X(3) *X(9)+X(4)*X(10) ) / (X(13) *(BC**2) ) 

F(14)=-(X(4)*X(9)-X(3)*X(10))/£S 

E=X(15) *BC  $  ES=E**2 

F(15)=(X(5)*X(11)+X(6)*X(12))/(X(15)*(BC**2)) 
F(16)=-(X(6)*X(11)-X(5)*X(12))/ES 
21  Dl=(l .-(X(13)**2) ) 

D2=X(13)/ (Dl**2) 

D3=( ( 16 . /SPI)-1 . ) 

D4=(D3*(X(13)**2)+1.) 

DO  46  1=13,14 
K  1+2 

F(I)=F(I) /PSDL 
46  F (K) =F (K) /PSDR 
C  BEARING  FORCE  CALCULATION 

FEL=~( ( (1 .-2.*F(14))*2.*(X(15)**2)/(D1**2) )+(PI*F(13)* 
l(l.+2.*(X(13)**2))/(Dl**2.5))) 
FPL=-((2.*F(14)-1.)*PI*X(13)/(2.*(D1**1.5))) 
l-(4.*F(13)*X(13)/(Dl**2))) 

FM=DVL*PSDL*BR/ (2 ,*(BC**2) ) 

FXL=( (FEL*SIN(X(14) )+FPL*COS (X(14) ) )*FM)-FXLS 
FYL=((-FEL*COS(X(14))+FPL*SIN(X(14)))*FM)-FYLS 
Dl=(l.-(X(15) **2) ) 

D2=X(15)/(D1**2) 

D3=((16./SPI)-1 .) 

D4=(D3*(X(15)**2)+1.) 

FER— (((l.-2.*F(16))*2.*(X(15)**2)/(Dl**2))+(PI*F(15)* 
1(1.+2.*(X(15)**2))/(D1**2.5))) 
FPR=-(((2.*F(16)-1.)*PI*X(15)/(2.*(D1**1.5))) 
l-(4.*F(15)*X(15)/ (Dl**2))) 

FM=DVR*PSDR*BR/ (2 . *(BC**2) ) 

FXR= ( ( FER*S IN (X ( 1 6) ) +FPR*COS (X ( 1 6 ) ) ) *FM) - FXRS 
FYR=((-FER*COS(X(16))+FPR*SIN(X(16)))*FM)-FYRS 
DO  49  1=13,14 
K=I+2 

F(I)=F(I)*PSDL 
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49  F(K)=F(K)*PSDR 
DO  57  1-1,6 

57  X(I)=X(I)-XS(I) 

STATE  VARIABLE  INCREMENTS 
F(1)=X(7) 

F(2)=X(8) 

F(3)=X(9) 

F(4)=X(10) 

F(5)-X(ll) 

F(6)-X(12) 

F(7)—(-SKX*(X(l)-X(3))-SKX*(X(l)-X(5))/DM 

F(8)=(-(X(2)-X(4) )-(X(2)-X(6) ))*SKY/DM 

F(9)«(-SKX*(X(3)-X(1))+FXL)/BM 

F(10)*(-SKY*(X(4)-X(2))+FYL)/BM 

F(U)-(-SKX*(X(5)-X(l))+FXR)/BM 

F(12)=(-SKY*(X(2))+FYR)/BM 

DO  58  1-1,6 

58  X(I)=X(I)+XS(I) 

RETURN 

END 

SUBROUTINE  RK4(F,Y,X,NN,HH) 

4TH  ORDER  RUNGE-KUTTA 
DIMENSION  Y(20),F(20) 

DIMENSION  SAVEY(20) , PHI (20) 

N=NN  $  H=HH 
M=0 

1  M=MH 

CALL  DERIV(F,Y,X) 

GO  TO  (2, 3, 4, 5) ,M 

2  DO  22  J=1,N 
SAVEY(J)=Y(J) 

PHI(J)-F(J) 

22  Y(J)-Y(J)+.5*F(J)*H 
X-X+. 5*H 
GO  TO  1 

3  DO  33  J-l.N 
PHI(J)=PHI(J)+2.0*F(J) 

33  Y(J)«SAVEY(J)+.5*H*F(J) 

GO  TO  1 

4  DO  44  J-l ,N 
PHI(J)=PHI(J)+2.0*F(J) 

44  Y (J)-SAVEY (J)+H*F (J) 

X-X+. 5*H 
GO  TO  1 

5  DO  55  J-l.N 

55  Y(J)-SAVEY(J)+(PHI(J)+F(J))*H/ 6.0 
RETURN 
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A. 2  Linear  Rotor  Bearing  System  Model  Programming 

PROGRAM  RWE4 ( RWE4 , OUTPUT , TAPE5'RWE4 , TAPE6=0UTPUT) 

DIMENSION  A(40,40) ,XR(40) ,XI(40) ,BRE(40,40) ,BI(40,40) 

C  THIS  PROGRAM  GENERATES  EIGENVALUES  AND  EIGENVECTORS  FOR  THE  MYRICK 
C  ROTOR  BEARING  SYSTEM  , 

C  SKX=SHAFT  STIFFNESS  IN  X  DIRECTION=5348 . 6  LBF/IN" 

C  SKY=SHAFT  STIFFNESS  IN  Y  DIRECTION' 5348 . 6  LBF/IN^ 

C  BM= JOURNAL  MASS=.0096  LBFSECZ/IN 

C  DM'DISC  MASS'. 0875  LBFSECZ/IN 

C  W  =  h( JOURNAL  +  DISC)  WEIGHT=20 . 625  LBF 

C  DVL=B EARING  FLUID  DYNAMIC  VISCOSITY  (LBFSEC/INZ) 

C  CL= JOURNAL  LENGTH=1.603  IN 

C  BR= JOURNAL  RADIUS' 1.081  IN 

C  BC'JOURNAL  TO  BEARING  RADIAL  CLEARANCE' . 003  IN 

READ (5, 11)SKX,SKY,BM,DM,W,CL,BR,BC 
PI'3. 1415927 
SPI=9. 8696044 
DVL=5.0E-06  $  DVR=5 .0E-06 
C  SET  INITIAL  ECCENTRICITY  RATIO 

ERR=. 075  $  ERL'. 075 
PRINT  1 1 ,DM,BM 
PRINT  11 

21  Dl'(l .-(ERR**2) ) 

D2=ERR/(D1**2) 

D3=( (16. /SPI)-1 . ) 

D4=(D3*(ERR**2)+1.) 

C  CALCULATE  ANGULAR  VELOCITY 

PSDR=-4.*W*(BC**2)*(D1**2)/(PI*DVR*CL*ERR*SQRT(D4)*BR) 
FM'DVR*PSDR*BR/ (BC**2) 

FMD'FM/PSDR 

PER=ATAN( (PI/4. )*SQRT(D1) /ERR) 

C  CALCULATE  BEARING  COEFFICIENTS 

BKEER'FM* (-2. *D2-(D2*(ERR**2) *4. /Dl) ) /BC 
BKEPR=-PI*FM*ERR/ (4 . * (Dl**l . 5) *ERR*BC) 

BDEER=-PI*FMD*(1 .+2 . * (ERR**2) )/(2.*(Dl**2. 5)*BC) 

BDEPR=2 . *ERR*D2*FMD/ (ERR*BC) 

BKPER'FM*((PI/(4.*(D1**1.5)))+(.75*PI*(ERR**2)/((D1**2.5)))) 

1/BC 

BKPPR=-ERR*D2*FM/ (ERR*BC) 

BDPER=2 . * (D2/BC) *FMD 

BDPPR=-PI*FMD*ERR/ (2 . *(D1**1 . 5)*ERR*BC) 

S-SIN(PER) 

C'COS (PER) 

C  BEARING  COEFFICIENTS  IN  X-Y  FORM 

BKXXR=((BKEER*S+BKEPR*C)*S+(BKPER*S+BKPPR*C)*C) 

BKXYR=((-BKEER*C+BKEPR*S)*S+(-BKPER*C+BKPPR*S)*C) 

BKYXR=(-(BKEER*S+BKEPR*C)*C+(BKPER*S+BKPPR*C)*S) 

BKYYR= ( - ( -BKEER*C+BKEPR*S ) *C+(-BKPER*C+BKPPR*S) *S ) 

BOXXR=((BDEER*S+BDEPR*C)*S+(BDPER*S+BDPPR*C)*C) 

BDXYR=((-BDEER*C+BDEPR*S)*S+(-BDPER*C+BDPPR*S)*C) 

BDYXR' ( - ( BDEER*S+BDEPR*C) *C+(BDPER*S+BDPPR*C) *S) 
BDYYR=(-(-BDEER*C+BDEPR*S)*C+(-BDPER*C+BDPPR*S)*S) 
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NC=0  $  N=12 

C  INITIALIZE  A  MATRIX 
DO  26  1=1, N 
DO  27  J=1,N 
27  A(I, J)=0.0 
26  CONTINUE 
DO  3  1=1,6 
KI=I+6 
A(I,KI)=1 . 

3  CONTINUE 

C  LOAD  A  MATRIX 

A(7 , 1)=-(SKX*2 . ) /DM 

A(7,3)=SKX/DM 

A(7,5)=SKX/DM 

A(8,2)=-(SKY*2.)/DM 

A(8 ,4)=SKY/DM 

A(8,6)=SKY/DM 

A(9,1)=SKX/BM 

A ( 9 , 3 ) =- ( SKX+BKXXR) /BM 

A(9,4)=-BKXYR/BM 

A(9 ,9)=-BDXXR/BM 

A ( 9 , 1 0) =-BDXYR/ BM 

A(I0, 2)=SKY/BM 

A(10,3)=-BKYXR/BM 

A ( 1 0 , 4 ) =- ( SKY+BKYYR) / BM 

A(10,9)=-BDXYR/BM 

A(10,10)=-BDYYR/BM 

A(ll , 1)=SKX/BM 

A ( 1 1 , 5) =-( SKX+BKXXR) / BM 

A(11,6)=-BKXYR/BM 

A(11,11)=-BDXXR/BM 

A(11,12)=-BDXYR/BM 

A(12,2)=SKY/BM 

A ( 1 2 , 5 ) =-BKYXR/ BM 

A(  1 2 , 6) =- (SKY+BKYYR) /BM 

A(12,11)=-BDYXR/BM 

A(12,12)=-BDYYR/BM 

DO  89  11=1,12 

DO  88  JJ=1 , 12 

88  A(II, JJ)*-A(II, JJ) 

89  CONTINUE 
J=1 

C  CALCULATE  EIGENVALUES  AND  EIGENVECTORS 

CALL  QRHMOD(A,N, J,XR,XI,BRE,BI) 
PSDR=-PSDR 
PRINT  2 , PSDR , ERR , PER 

C  PRINT  EIGENVALUES 

PRINT  30, (XR(I) ,XI(I) , 1=1 ,N) 

PRINT  11 
PRINT  11 

C  INCREMENT  ECCENTRICITY  RATIO 

ERR=ERR- .0001 
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A.  3 
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C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


Linear  Rotor  Bearing  System  and  Labyrinth  Seal  -  Coupled  (RWE4S) 
and  Isolated  Seal  and  Disc  (RWE4P)  Programming 

PROGRAM  RWE4S (RWE4S .OUTPUT .TAPE5-RWE4S .TAPE6-0UTPUT) 

DIMENSION  AM(40,40),XR(40) ,BRE(40,40) ,BI (40 ,40) ,PC(4,4) 

DATA  PC/ 16*0./ 

THIS  PROGRAM  GENERATES  EIGENVALUES  AND  EIGENVECTORS  FOR  A  ROTOR 
BEARING  SYSTEM  WITH  SINGLE  LABYRINTH  SEAL 
SKX-  X  DIRECTION  ROTOR  ELASTICITY-5348.6  LBF/IN 
SKY-  Y  DIRECTION  ROTOR  ELASTICITY-5348.6  LBF/IN 
BM  -  JOURNAL  MASS-. 0096  LBF*SEC2/IN 
DM  -  DISC  MASS-. 0879  LBF*SEC2/IN 

W  -  h  WEIGHT  OF  ROTOR  BEARING  SYSTEM  ROTATING  PARTS-20.625  LBF 

CL  -  JOURNAL  LENGTH- 1.000  IN 

BR  -  JOURNAL  RADIUS- 1.081  IN 

BC  -  RADIAL  BEARING  CLEARANCE-. 003  IN 

CLC-  LABYRINTH  SEAL  AXIAL  LENGTH 

WD  -  LABYRINTH  SEAL  CIRCUMFERENTIAL  LENGTH-PI*RD*2. 

PSI-  UPSTREAM  SEAL  PRESSURE 
TI  -  SEAL  TEMPERATURE- 5 6 OR 
R  -  GAS  CONSTANT  (AIR)-247104  IN2/SEC2R 
A  =  UPSTREAM  SEAL  STRIP  HEIGHT 
B  -  DOWNSTREAM  SEAL  STRIP  HEIGHT 
DC  -  DISC  TO  SEAL  RADIAL  CLEARANCE 
AMNI-  INITIAL  UPSTREAM  MACH#  (ASSUMED) 

CK  -  SPECIFIC  HEAT  RATIO 

CWB-  DOWNSTREAM  SEAL  STRIP  DISCHARGE  COEFFICIENT 
PSE-  SEAL  DISCHARGE  PRESSURE 
DV  -  SEAL  FLUID  DYNAMIC  VISCOSITY 
SPI-  PI**2 

ERR, ERL-  JOURNAL  ECCENTRICITY  RATIOS,  RIGHT&LEFT  JOURNALS 

RD  -  DISC  RADIUS 

VK  -  KINEMATIC  VISCOSITY 

READ (5 ,11) SKX,SKY,BM,DM,W,CL,BR,BC,CLC,WD,PSI,TI,R,A,B 
1 , DC , AMNI , CK , RD , CWB , PSE 
DV-2.8E  09 
PI-3.1415927 
SPI-9. 8696044 

CK1- (CK- 1 . ) /CK  $  CK1R-1./CK1  $  CKR-l./CK 
CK2- (CK+1. ) / (2.*(CK-1. ) ) 

CK3- (3. -CK) / (2. *(CK- 1 . ) ) 

CK4-(CK-1. )/2.  $  CK5-(CK+l.)/2.  $  CK5R-1./CK5 
CK6-CK2*2.  $  CK7-(1.-2.*CK)/CK 
CK8- SQRT ( 1 . +CK4 ) 

CK9-(1 . / CK5) **CK1R 
CK10- 1 . / ( ( 1 . +CK4 ) **CK2 ) 

CK1 1-SQRT( (CK/R) * (CK5R**CK6) ) 

BEARING  FLUID  DYNAMIC  VISCOSITY 
DVL-5 . 0E-06  $  DVR-5.0E-06 

INITIALIZED  JOURNAL  ECCENTRICITY  RATIO  (VALUE  AS  REQUIRED) 

ERR-. 076 

ITERATE  SEAL  PARAMETERS 
DC  22  IBI-1,10 
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C  SEAL  ENTRANCE/EXIT  AREAS 

21  AB=(DC-B)*WD 
AA-(DC-A)*WD 

C  CALCULATE  SEAL  CHAMBER  PRESSURE 

AMNS=AMNI 

24  FMAS=AMNS/((1.+CK4*(AMNS**2))**CK2) 

PNP-1 .  /  (  ( 1 .+CK4* (AMNS**2) ) **CK1R) 

CWA*-.3179464*PNP+. 9123729 
IF(AMNS . LE. 0 . )AMNS=0 . 

IF(AMNS.LE.O.)GO  TO  37 
PN-PNP*PSI 

AMN=SQRT((1./CK4)*(((PSI/PN)**CK1)-1.)) 

CWB  =*-.31 79464*  (PSE/PN)+.  9 123729 
BMNS-SQRT ( ( 1 . /CK4) * ( ( (PN/PSE) **CK1 ) -1 . ) ) 

IF(BMNS.GE.1)BMNS=1. 

IF (BMNS . GE1 . ) CWB- . 74 

FMBS-BMNS/ ( (1 .+CK4* (BMNS**2) ) **CK2) 

PNN°AA*PSI*FMAS*CWA/ (AB*FMBS*CWB) 

IF (BMNS . GE . 1 . ) PNN-PSE/ CK9 
IF (PNN . GT . PN) AMNS-AMNS- . 00001 
IF (PNN . LT . PN) AMNS-AMNS+. 0000 1 

IF(PNN.LE.(PN+.0005*PN) .AND. PNN. GE. (PN-.0005*PN)  GO  TO  37 
IF(AMN.GE.  1 .  )AMNS**1 . 

GO  TO  24 

37  IK-1  $  N-16  $  NN-20  $  MM-10  $  NP-4 

PRINT  11,AMN,AMNS,BMNS,PN,PNN 
PRINT  11 
AMNI-AMN 
D1=(1.-(ERR**2)) 

D2-ERR/ (Dl**2) 

D3-((16./SPI)-1.) 

D4-(D3*(ERR**2)+1.) 

C  CALCULATE  ROTOR  ANGULAR  VELOCITY 

PSDR--4 . *W* (BC**2) * (Dl**2) / (PI*DVR*CL*ERR*SQRT(D4)*BR) 
FM=DVR*PSDR*BR/ (BC**2) 

FMD-FM/PSDR 

C  CALCULATE  JOURNAL  ATTITUDE  ANGLE 

PER-ATAN( (PI/4 . ) *SQRT (Dl) /ERR) 

C  CALCULATE  BEARING  COEFFICIENTS 

BKEER»FM*(-2 . *D2-(D2*(ERR**2) *4 . /Dl ) ) /BC 
BKEPR— PI*FM*ERR/ (4 .*(D1**1 . 5)*ERR*BC) 

BDEER— PI*FMD* (l.+2.*(ERR**2))/(2.*(Dl**2.5)*BC) 

BDEPR-2 . *ERR*D2*FMD/ ( ERR*BC) 

BKPER-FM*((PI/(4.*(Dl**1.5)))+(.75*PI*(ERR**2)/((Dl**2.5)))) 

1/BC 

BKPPR— ERR*D2*FM/ (ERR*BC) 

BDPER-2 .*(D2/BC)*FMD 

BDPPR— PI*FMD*ERR/ ( 2 . * ( Dl ** 1 . 5) *ERR*BC) 

S-SIN(PER) 

C-COS(PER) 

BKXXR-((BKEER*S+BKEPR*C)*S+(BKPER*S+BKPPR*C)*C) 

BKXYR-((-BKEER*C+BKEPR*S)*S+(-BKPER*C+BKPPR*S)*C) 

BKYXR- (  - ( BKEER*S+BKEPR*C ) *C+(BKPER*S+BKPPR*C ) *  S ) 


! 
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BKYYR-(  - ( -BKEER*C+BKEPR* S ) *C+ ( -BKPER* C+BKPPR* S ) * S ) 

BDXXR-((BDEER*S+BDEPR*C)*S+(BDPER*S+BDPPR*C)*C) 

BDXYR»(-BDEER*C+BDEPR*S)*S+(-BDPER*C+BDPPR*S)*C) 

BDYXR*(  -(BDEER*S+BDEPR*C)*C+(BDPER*S+BDPPR*C)*S) 

BDYYR- (  - ( -BD EER* C+BD EPR* S ) *C+( -BDPER*C+BDPPR*S ) *S) 

ACH-CLC*WD 

ED-ERR*BC 

C  ZERO  A  MATRIX 
DO  26  I-l.N 
DO  27  J-1,N 
27  AM(I, J)-0.0 
26  CONTINUE 

C  XI  THRU  X6  INITIALIZATION  IN  AM 
DO  3  1-1,6 
KI-I+6 
AM(I,KI)-1. 

3  CONTINUE 
RD2-RD/2. 

PSDU— PSDR*RD 
HN-DC 

AA-(HN-A) *WD 
AB- (HN-B) *WD 

C  SEAL  FRICTION  COEFFICIENTS 

UN=(PSDU/2.) 

C  REO-REYNOLDS#  AT  SEAL  BASE 

REO=ABS (PN*UN*HN/ (2 . *R*TI*DV) ) 

FFEO=-l./4.38 
AFO=. 04176 

IF(REO.GT.l.E+05)FFEO=-l./7.4 
IF(RE0.GT.1.E+05)AF0=. 01423 
C  REI=REYNOLDS#  AT  DISC 

REI-ABS (PN* (PSDU-UN) *HN/ ( 2 . *R*TI*DV) ) 

FFEI--1./4.48 

AFI-.04176 

IF(REI . GT . 1 . E+05 ) FFEI=- 1 . / 7 . 4 
IF(REI.GT.1.E+05)AFI». 01423 
FFNI— l./FFEI 
FFNO-l./FFEO 

C  DISC  AND  SEAL  BASE  FRICTION  FACTORS 
FFI-(REI**FFEI)*AFI 
FFO-(REO**FFEO) *AFO 
RT-R*TI 

C  SEAL  AXIAL  FLOW  COEFFICIENTS 
BK-SQRT (CK/ RT) *PSE 
C  DISCHARGE  COEFFICIENT  AT  B 

CWBN—.3179464*(PSE/PN)+.  9123729 
CWBK- . 31 79464*PSE/ (PN**2) 

C  MACH^  . LE. 1 

IF( (PSE/PN) . LE. CK9) CWBN- . 74 
IF( (PSE/PN) . LE . CK9) CWBK=0 . 

C  MACH#  AT  B  .LE.l 

BMN-SQRT( ( 1 . /CK4) *( ( (PN/PSE)**CK1 )-l . ) ) 

IF( (PSE/PN) .LE.CK9)BMN-1. 
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GMB 1- ( 1 .+CK4* (BMN**  2 ) ) 

GMB-SQRT(GMBl) 

FMB«BMN*GMB 

FMBK- ( (CK4* ( (BMN**2) )+GMBl) /GMB) * (PSE**CK7) / 

1 ( (PN**CKR) *CK*BMN) 

C  AXIAL  FLOW  AT  B-COEFFICIENTS 

PMBK-BK*AB* (CWBN*FMBK+FMB*CWBK) 

HMBK»BK*CWBN*FMB*WD 
C  DISCHARGE  COEFFICIENT  AT  A 

CWAN*-. 3 179464* (PN/PSI) +. 9123729 
AK-SQRT (CK/ (R*TI) )*PSI 
CWAK—  .3179464/PSI 
IF( (PN/PSI) . LE.CK9)CWAN*. 74 
IF( (PN/PSI) .LE.CK9)CWAK-0. 

AMN*SQRT( ( 1 . /CK4) * ( ( (PSI/PN) **CK1 ) -1 . ) ) 

IF( (PN/PSI) .LE.CK9)AMN-1. 

GMA1-(1 .+CK4* (AMN**2) ) 

GMA-GMA1**CK2 

FMA-AMN/GMA 

FMAK“-(GMA-(AMN**2) *CK5*(GMA1**CK3) ) * (PN**CK7) / ( (GMA1**CK6) 

1* (PSI**CKR) *CK*AMN) 

C  AXIAL  FLOW  AT  A-COEFFICIENTS 

PMAK-AK*AA* ( CWAN*  FMAK+FMA*CWAK) 

HMAK“CWAN*FMA*WD*AK 
C  CONTINUITY  EQUATION  COEFFICIENTS 

C1PT-CLC*HN 
CIP-(PMAK-PMBK) *RT/WD 
C1H- (HMAK-HMBK) *RT/WD 
C1US— CLC*PN*HN 
C1PS-CLC*HN*UN 
C1HT— CLC*PN 
C1HS— CLC*PN*UN 
FFOE* (FFNO+1 . ) *FFEO 
FFIE-(FFNI+1.)*FFEI 
VK«DV*RT/PN 

C  MOMENTUM  EQUATION  FRICTION  TERMS 

FIPK— (AFI/FFNI) * (REI**FFIE) * ( (PSDU-UN) **3) *HN*PN*CLC 
1/ (4 .*DV*(RT**2) ) 

FOPK-(AFO/FFNO)*(REO**FFOE) *(UN**3) *HN*PN* (CLC+2 . *HN) 

1/ (4 .*DV*(RT**2) ) 

FIHK“( (AFI*UN/FFNI) * ( ( (PSDU-UN) *HN/ ( 2 . *VK) ) **FFIE) ) 

1*PN* ( (PSDU-UN) **2) *CLC/ (4 . *RT*VK) 

FIUK- ( (AFI*HN/FFNI) * ( ( (PSDU-UN) *HN/ ( 2 . *VK) ) **FFIE) ) 

1*PN* ( (PSDU-UN) **2) *CLC/ ( 4 . *RT*VK) 

FOHK-( (AFO*UN/FFNO) *( (UN*HN/ (2 . *VK) ) **FFOE) ) 

1*PN*(UN**2) * (CLC+2 . *HN) / (4 . *RT*VK) 

FOUR- ( (AFO*HN/FFNO) * ( (UN*HN/ ( 2 . *VK) ) **FFOE) ) 

1*PN* (UN**2) * (CLC+2 . *HN) / (4 . *RT*VK) 

C  MOMENTUM  EQUATION  COEFFICIENTS 

C2P-(.5*(FFI*((PSDU-UN)**2)*CLC  FFO*(UN**2)*(CLC+2)*HN))/RT 
D+FIPK+FOPK 

C2H“(-FFO*FN*(UN**2)/RT)+FIHK+FOHK 

C2U“(-FFI*PN*( PSDU-UN) *CLC/RT)+(-FFO*PN*UN* (CLC+2 . *HN) /RT)+FIUK+ 
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1F0UK 

C2PT— CLC*UN*HN/RT 
C2PS«-CLC*HN*  (  (  (UN**2)  /RT-1 . ) 

C  2UT* -CLC*PN*HN/ RT 
C2US— CLC*PN*2 .  *UN*HN/RT 
C  C2HT-- CLC*PN*UN/RT 

C2HS— CLC*PN*  ( 1 . +(UN**2)  /RT 
C  ZERO  SEAL  MATRIX 

DO  100  IP-l.NP 
DO  101  JP-l.NP 
101  PC(IP,JP)»0 
100  CONTINUE 

C  SEAL  MATRIX  LOADING 

PC(1,1)“C1P/C1PT 
PC(1,2)«C1PS/(RD*C1PT) 

PC(1,4)-C1US/(RD*C1PT) 

PC(2,1)=!-C1PS/(RD*C1PT) 

PC(2, 2)*C1P/C1PT 
PC(2,3)=-C1US/(RD*C1PT) 

PC(3,1)  — (C2P+(C1P*C2PT/C1PT))/C2UT 
PC(3,2)=-(C2PS+(C1PS*C2PT/C1PT))/(RD*C2UT) 

PC ( 3 , 3) ”-C2U/ C2UT 

PC(3,4)— (C2US+(C1US*C2PT/C1PT))/(RD*C2UT) 

PC(4 , 1) =(C2PS+(C1PS*C2PT/C1PT) ) / (RD*C2UT) 

PC(4 , 2) — (C2P+(C1P*C2PT/C1PT) )/C2UT 

PC(4 , 3) =(C2US+(C1US*C2PT/C1PT) ) / (RD*C2UT) 

PC(4,4)--C2U/C2UT 

C  DISC  TO  SEAL  COUPLING  COEFFICIENTS 
CAH— C1H/C1PT 

CCH»(C1H*C2PT/ (C1PT*C2UT) )+C2H/C2UT 
CBH-C1HS/(RD*C1PT) 

CDH— (C2PT*C1HS/ (RD*C1PT*C2UT) )-(C2HS/ (RD*C2UT) ) 
CAHD— C1HT/C1PT 

CCHD-(C2PT*C1HT/ (C1PT*C2UT) )+(C2HT/C2UT) 

KP-13 

KP3-KP+3 

C  A  MATRIX  LOADING  WITH  SEAL  COEFFICIENTS 
DO  50  I-KP.KP3 
DO  51  J-KP.KP3 
KPJ-J-KP+1 
KPI-I-KP+1 

51  AM(I,J)=PC(KPI,KPJ) 

50  CONTINUE 

C  SEAL  TO  DISC  COUPLING 
AM(7 , 14) -PI*RD*CLC/DM 
AM( 8 , 1 3) »PI*RD*CLC/ DM 
C  DISC  TO  SEAL  COUPLING 
AM(13,2)-CAH 
AM(15,2)-CCH 
AM(14,2)-CBH 
AM(16,2)-CDH 
AM(13,8)-CAHD 
AM(15,8)-CCHD 


*4  <**.,-#* 
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C  A  MATRIX  LOADING  WITH  ROTOR  BEARING  SYSTEM  COEFFICIENTS 
AM(7,l)-(-(SKX*2.)/DM) 

AM(7,3)«SKX/DM 

AM(7,5)-SKX/DM 

AM(8,2)»(-(SKY*2.)/DM) 

AM(8,4)-SKY/DM 
AM(8,6)«SKY/DM 
AM(9,1)-SKX/BM 
AM(  9, 3)  — (SKX+BKXXR) /BM 
AM(9,4)— BKXYR/BM 
AM(9, 9)— BDXXR/BM 
AM(9,10)— BDXYR/BM 
AM(10,2)»SKY/BM 
AM(10,3)  —  BKYXR/BM 
AM( 1 0 , 4) — ( SKY+BKYYR) /BM 
AM(10,9)  —  BDYXR/BM 
AM(10,10)—  BDYYR/BM 
AM(11,1)-SKX/BM 
AM(1 1,5)— (SKX+BKXXR)  /BM 
AM(11,6)—  BKXYR/BM 
AM(  1 1 , 1 1) —BDXXR/BM 
AM(  11,12) — BDXYR/BM 
AM(12,2)-SKY/BM 
AM(12, 5)— BKYXR/BM 
AM(  1 2 , 6  )—(  SKY+BKYYR) /BM 
AM(  12, 11)—  BDYXR/BM 
AM(12, 12)—  BDYYR/BM 
73  DO  89  11-1,12 
DO  88  JJ-l.N 

88  AM(II.JJ)—  AM(II.JJ) 

89  CONTINUE 
J-l 

C  CALCULATE  EIGENVALUES  &  EIGENVECTORS  FOR  AM 
C  RWE4P  ONLY 

CALL  QRHMOD (PC,NP,J,XR,XI,BRE,BI) 

C  RWE4S  ONLY 

CALL  QRHMOD (AM, N, J,XR, XI, BRE.BI) 

PSD— PSDR 

PRINT  11,A,B,PSI,PSE 
PRINT  6 , CLC , RD , DC , REI , REO , FFI , FFO 
PRINT  2, PSD, ERR, PER 
C  PRINT  EIGENVALUES 

PRINT  30, (XR(I) ,XI(I) ,1-1 ,N) 

PRINT  11 

C  SEAL  PARAMETERS  VARIATIONS 

A-A+.0005 
CLC -CLC+. 0249 
RD-RD+1. 

PSI-PSI+. 1549 
DC-DC+.0802 
22  CONTINUE 
C  GO  TO  23 
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C 


C 


SPEED  VARIATIONS 
ERR- ERR-. 01 

IF(ERR.LT. .074)GO  TO  23 
IF(ERR.GT. .5) GO  TO  23 


GO  TO  21 


23  CONTINUE 

PRINT  EIGENVECTORS 

PRINT  30,((B*B(I.J  •JJ5VI:1.™  1 

2  FORMAT (9X,'PSDR  -  .E12.4.3X,  ERR 

in  FORMAT(5X,E15.6,5X,E15.6) 


11  FORMAT (El 5. 6) 

6  FORMAT(10X,7E15.5/) 


J-l.N) 

:  ’ ,E12.4,3X, 'PER  - 


STOP 

END 


' ,E12.4/) 


APPENDIX  B 


B.l  Beaman  13^  Bearing  Coefficients  (Equation  3.25) 
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B.2  XY  Bearing  Coefficients  (Equation  3.28) 
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(C«e  Cos  4  4"  ^  cos 

(  Kst.  c-o$  4  +■  ^ 
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B.3  Incremental  Seal  Strip  Mass  Flow  Rate  Coefficients 
From  Equation  3.49e: 
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